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ABSTRACT 


Numerical  weather  prediction  models  have  not  produced  accurate 
precipitation  forecasts,  especially  short-term  forecasts  of  significant 
precipitation  events.  One  reason  for  this  has  been  that  numerical 
models  are  normally  initialized  with  nondivergent  winds.  This  means 
the  model  must  develop  a  vertical  motion  field  and  an  associated 
precipitation  field.  Therefore,  the  initial  precipitation  rate  is 
underforecast  and  the  precipitation  forecast  itself  is  adversely 
affected.  One  method  of  solving  this  problem  is  the  initialization 
of  the  divergent  wind  component  (hereafter  termed  divergent 
initialization).  Previous  divergent  initialization  attempts  have  been 
primarily  on  the  synoptic  scale.  These  attempts  did  not  produce  a 
significant  change  in  the  initial  precipitation  rate.  Divergent 
initialization  on  the  mesoscale  will  be  attempted  here. 

The  general  divergent  initialization  procedure  proceeds  as 
follows:  omega  values  are  diagnosed  using  the  omega  equation; 
velocity  potentials  are  derived  from  the  vertical  velocities  with  the 
continuity  equation;  the  divergent  wind  components  are  obtained  from 
the  velocity  potentials;  geopotentials  are  calculated  on  sigma 
surfaces  using  a  balance  equation  with  contributions  from  both  the 
nondivergent  and  divergent  wind  components;  finally,  balanced 
temperatures  are  derived  using  the  hydrostatic  equation. 

A  scale  analysis  was  performed  on  the  vertical  velocity  and 
divergence  equations  to  determine  the  forms  appropriate  for  the 
mesoscale  (grid  increments  from  60  to  200  km).  The  scale-dependent 

i 


iv 


differences  of  divergent  initialization  on  the  synoptic  and  mesoscales 
were  studied.  Both  nondivergent  and  divergent  wind  components  are 
required  for  these  equations.  To  obtain  them,  boundary  conditions  on 
stream  function  and  velocity  potential  are  required.  The  low-wave- 
number  boundary  variation  of  stream  function  and  velocity  potential 
should  be  specified  accurately  to  minimize  the  influence  of  the 
boundary  values  on  the  solution  in  the  domain  interior.  A  method  for 
specifying  boundary  conditions  on  velocity  potential  with  accurate 
low-wavenumber  boundary  variation  was  developed. 

The  forecast  model  used  was  developed  at  The  Pennsylvania  State 
University.  The  version  used  here  had  six  levels,  low- resolution 
planetary  boundary  layer  physics,  and  a  grid  increment  of  120  km. 

The  mesoscale  omega  equation  was  solved  by  three-dimensional 
relaxation  over  the  domain.  The  observed  rainfall  rate  was  used  to 
construct  a  parabolic  omega  profile.  A  heating  rate  was  derived  from 
the  profile  and  used  as  input  for  the  diabatic  term  in  the  omega 
equation.  On  the  mesoscale,  the  largest  single  term  in  the  omega 
equation  was  the  diabatic  term.  The  greatest  uncertainty  in  the 
calculation  of  omega  values  was  the  representativeness  of  the  observed 
precipitation  rate. 

Five  12-hour  forecasts  were  conducted,  two  with  unbalanced  initial 
conditions  and  three  with  different  balanced  initial  conditions. 

Three  forecasts  had  essentially  nondivergent  initial  conditions.  The 
fifth  forecast  was  balanced  on  sigma  surfaces  and  was  initialized  with 


the  total  wind,  the  divergent  part  of  which  was  derived  from  the 
diagnosed  omega  field.  The  noise  characteristics  of  the  five 
forecasts  did  not  differ  greatly  from  those  reported  by  other 
investigators.  However,  the  forecast  from  the  divergent  initialization 
produced  more  precipitation  than  the  other  forecasts,  especially  in 
the  first  three  hours.  Also,  the  divergent  initialization  produced 
a  more  evenly  distributed  precipitation  prediction.  The  forecast 
model  retained  the  initial  divergence.  The  effect  of  divergent 
initialization  was  significant  for  this  case  and  would  likely  improve 
short-range  precipitation  forecasts  from  mesoscale  numerical  weather 


prediction  models. 
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1.0  INTRODUCTION 


Perhaps  the  four  most  important  events  to  influence  meteorology 
in  the  twentieth  century  are: 

(1)  the  development  of  the  concept  of  a  front, 

(2)  the  development  of  a  global  rawinsonde  network  for  the 
twice-daily  sampling  of  the  atmosphere, 

(3)  the  invention  of  the  electronic  computer,  and 

(4)  the  use  of  earth-orbiting  satellite  platforms  for  the 
observation  of  the  atmosphere. 

Simultaneously,  meteorologists  have  gained  increased  knowledge  and 
insight  into  the  physics  of  the  atmosphere.  All  of  this,  coupled 
with  the  very  rapid  developments  in  computer  technology,  has  led  to  an 
ever  increasing  interest  in  the  mathematical  modeling  of  atmospheric 
structures  in  general  and  in  numerical  weather  prediction  (NWP)  in 
particular. 

Primitive  equation  models  of  the  atmosphere  require  initial 
conditions  for  the  mass  and  momentum  variables.  The  initial  conditions 
are  derived  from  a  limited  number  of  observations,  particularly  at 
synoptic  times.  Many  observations  which  are  available,  including  some 
at  synoptic  times,  are  not  used  to  assist  in  defining  the  initial 
conditions.  Additionally,  most  observations  that  are  used  contain 
errors.  After  a  meteorological  field  is  analyzed,  it  must  then  be 
interpolated  vertically  and/or  horizontally  to  define  the  field  at 
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the  model  grid  points.  This  interpolation  process  also  introduces 
errors.  Therefore,  the  entire  data  analysis  process  introduces  errors 
even  though  the  process  itself  is  ersential  for  a  model  integration. 

Once  the  model  variables  have  been  analyzed  to  the  model  grid 
points,  it  is  usually  desirable  to  balance  the  mass  and  momentum  fields 
Balancing  means  that  the  mass  and  momentum  fields  are  made  compatible 
with  one  another  so  that  the  nonmeteorological  accelerations  produced 
at  the  beginning  of  a  forecast  are  minimized.  Normally,  balancing  is 
accomplished  by  requiring  that  the  mass  and  momentum  fields  initially 
obey  some  diagnostic  relationship.  Examples  of  diagnostic  relationship 
include  the  geostrophic  wind  equation,  the  gradient  wind  equation,  and 
the  balance  equation.  An  initialization  procedure  of  this  kind,  i.e., 
one  which  uses  data  primarily  from  one  time  and  that  does  not  use  the 
model  equations  for  balancing,  is  known  as  a  static  initialization  (SI) 

For  complex  baroclinic  models  that  contain  complex  parameteriza- 
tlons  of  diabatic,  frictional,  and  radiational  effects,  it  is  usually 
difficult  if  not  impossible  to  find  simple  but  adequate  diagnostic 
relationships  between  the  mass  and  momentum  fields.  This  fact  has 
led  to  a  new  kind  of  initialization  procedure  known  as  dynamic 
initialization  (DI).  The  distinguishing  feature  of  DI  is  that  the 
model  equations  themselves  are  used  for  the  balancing  step. 

In  a  standard  SI,  the  initial  conditions  are  customarily  non- 
divergent  because  of  the  diagnostic  balance  relationship  usually 
employed.  Although  they  are  relatively  simple,  nondivergent  initial 
conditions  do  not  permit  a  model  to  initially  forecast  precipitation 
at  the  observed  rate.  Since  mesoscale  NWP  models  will  be  used 
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primarily  for  short-range  (6-  to  24-hour)  forecasts,  it  is  especially 
desirable  for  the  model  to  accurately  forecast  the  initial  precipitation 
rate.  If  the  initial  precipitation  rate  is  greatly  underforecast , 
then  the  model  probably  will  not  forecast  precipitation  correctly  for 
periods  on  the  order  of  twelve  hours  because  the  model  should  not 
overforecast  precipitation  late  in  the  12-hour  period. 

A  method  of  correcting  this  precipitation  rate  problem,  previously 
employed  primarily  on  the  synoptic  scale,  has  been  the  addition  of  a 
divergent  component  to  the  nondivergent  initial  conditions.  Normally, 
this  procedure  consists  of  using  the  omega  equation  to  get  vertical 
velocities.  Then  the  continuity  equation  with  appropriate  boundary 
conditions  is  solved  to  get  the  velocity  potential  and  the  divergent 
wind  field.  The  divergent  wind  component  is  then  added  to  the 
nondivergent  initial  conditions.  Generally,  results  with  this  technique 
on  the  synoptic  scale  have  produced  little  change  over  results  based  on 
nondivergent  initial  conditions.  The  objective  of  this  thesis  is  to 
study  how  the  divergent  wind  component  can  be  incorporated  in  a  static 
initialization  on  the  mesoscale  in  such  a  way  as  to  improve  the 
initial  forecast  precipitation  rate.  This  should  in  turn  improve  the 
model's  forecast  precipitation  amounts  for  varying  time  periods. 

1.1  Review  of  the  numerical  model  initialization  problem 

Development  of  initialization  techniques  has  paralleled  the 
development  of  numerical  models  themselves.  Therefore,  it  is 


appropriate  to  review  briefly  the  history  of  NWP  before  discussing 
initialization. 

1.1.1  Brief  history  of  numerical  weather  prediction 

The  discovery  leading  to  modern  NWP  was  the  realization  that  the 
forecast  problem  can  be  treated  mathematically  as  an  initial  value¬ 
boundary  value  problem  for  the  nonlinear  partial  differential 
equations  which  govern  atmospheric  motion.  L.  F.  Richardson  (1922) 
was  the  first  to  apply  this  idea.  Richardson  encountered  several 
problems  which  we  are  careful  to  avoid  today.  One  of  the  problems 
that  led  to  his  failure  was  the  use  of  observed  rather  than  balanced 
data.  The  resultant  large-amplitude  nonmeteoro logical  waves  destroyed 
the  meteorological  forecast.  Additionally,  since  computers  were  not 
yet  invented,  a  large  amount  of  time  was  required  for  the  calculations. 
Richardson's  failure  resulted  in  a  lack  of  interest  in  NWP  that 
lasted  about  25  years. 

Renewed  interest  in  NWP  occurred  in  the  late  1940' s.  By  then 
the  use  of  rawinsondes  made  more  accurate  and  more  spatially  and 
temporally  numerous  observations  available.  Also,  the  electronic 
computer  had  been  invented.  Simple  "filtered"  models,  which  eliminated 
the  large-amplitude  fast  moving  waves  (and  thus  permitted  larger  time 
steps)  that  destroyed  Richardson's  forecast,  were  developed. 

As  computer  size  and  speed  increased,  primitive  equation  (PE) 
models  came  into  use  because  they  permit  greater  flexibility.  Gravity- 
inertia  waves,  which  are  meteorologically  important  in  the  atmosphere, 


5 


are  allowed  to  exist  in  PE  models.  The  philosophy  is  that  the  model 
which  most  closely  resembles  the  atmosphere  will  best  forecast  the 
atmosphere.  However,  this  added  flexibility  requires  much  greater 
care  in  the  formulation  of  the  initial  conditions.  One  aspect  of 
that  formulation  is  addressed  in  this  thesis. 

Improved  numerical  forecasts  are  hampered  by  three  constraints 
(Haltiner  and  Williams,  1975): 

(1)  The  temporally  and  spatially  continuous  atmosphere 
must  be  represented  by  data  at  discrete  points. 

Therefore,  the  entire  spectrum  of  atmospheric  motions 
cannot  be  resolved.  In  addition,  numerous  problems 
arise  due  to  the  approximation  of  derivatives  by 
finite  differences. 

(2)  We  do  not  completely  understand  the  physics  of  the 
atmosphere.  Processes  which  make  important  contributions 
to  the  evolution  of  the  atmosphere  are  not  adequately 
represented.  Also,  computers  are  limited  In  terms 

of  their  size  and  speed  in  spite  of  the  very  rapid 
advances  in  computer  technology. 

(3)  Inappropriate  specification  of  the  initial  conditions 
will  limit  forecast  accuracy  because  a  numerical 
forecast  is  an  initial  value  problem. 

Dutton  (1976)  pointed  out  that  even  if  the  initial  conditions 
are  perfect,  the  model's  ability  is  limited  because  imperfectly 
handled  nonlinear  interactions  lead  to  forecast  error  growth. 
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Nonlinear  interactions  produce  unresolvable  wavelengths  which 
erroneously  appear  on  the  grid  at  resolvable  scales  (this  phenomenon 
is  called  aliasing).  Because  the  nonlinear  interactions  are  always 
there,  the  error  increases  with  each  time  step.  Lorenz  (1965)  stated 
that  the  difference  between  the  forecast  and  actuality  will  eventually 
be  that  of  two  random  states  for  the  appropriate  time  of  day  and  year. 

1.1.2  Stages  of  the  initialization  process 

The  initialization  problem  has  defied  easy  solution.  The  basic 
stages  of  the  process  for  grid-point  models  are: 

(1)  observation, 

(2)  analysis  to  grid  points,  and 

(3)  balancing  the  mass  and  momentum  fields. 

The  first  stage  has  often  been  a  problem  because  the  observation 
density  is  usually  low.  The  observations  also  contain  human  and 
instrumental  errors  as  well  as  information  on  unresolvable  scales. 

The  second  stage  consists  of  interpolation  from  observation 
locations  to  grid-point  locations  and  to  a  fixed  time.  Cressman  (1959) 
stated  that  interpolation  should 

(1)  ensure  internal  consistency  among  variables, 

(2)  reduce  observation  errors,  and 

(3)  smooth  subgrid-scale  features. 
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The  third  stage  in  the  initialization  process  can  be  crucial. 
Thompson  (1961)  noted  that  large-scale  accelerations  are  due  primarily 
to  small  differences  between  forces.  That  is,  the  net  accelerations 
are  10  percent  of  the  magnitude  of  the  individual  accelerations. 
Therefore,  to  determine  the  accelerations  within  10  percent,  the  wind 
and  horizontal  pressure  gradient  force  must  be  known  within  1  percent. 
This  accuracy  has  not  been  attainable  with  most  present-day  observations 
so  observed  data  cannot  be  used  directly  as  initial  conditions. 

A  dynamic  balance  may  be  defined  as  follows:  mass  (temperature 
and  pressure)  and  momentum  fields  are  compatible  (in  "balance")  when 
the  only  accelerations  produced  by  the  mass  and  momentum  fields 
generate  features  of  meteorological  importance.  When  an  imbalance 
exists  between  the  mass  and  momentum  fields,  the  atmosphere  tends  to 
reduce  the  incompatibility  through  a  mutual  adjustment  process.  The 
theory  of  this  process  has  been  developed  by  Rossby  (1938),  Cahn  (1945), 
and  Washington  (1964),  among  others.  This  theory,  called  the  geostrophic 
adjustment  theory,  predicts  that  on  the  mesoscale,  the  wind  field 
dominates  the  height  field  in  the  adjustment  process. 

One  of  the  mechanisms  which  accomplishes  the  adjustment  process 
is  the  gravity-inertia  or  gravity  wave.  A  gravity  wave  is  a  wave 
driven  by  the  force  of  gravity  which  can  redistribute  mass  and  kinetic 
energy.  Gravity  waves  that  are  generated  by  imbalances  in  the  initial 
state  traverse  and  then  leave  the  domain  and  hence  usually  do  not  have 
a  major  influence  on  the  forecast  after  12  to  24  hours.  However, 
unrealistic  pressure  fluctuations  and  accompanying  vertical  motions 


t 


can  completely  obscure  meteorologically  important  motions  during  the 
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first  twelve  hours  (Sasaki,  1969;  Haltiner  and  Williams,  1975).  This 
fact  has  made  short-term  forecasting  difficult.  If  the  model  initial 
conditions  have  not  been  balanced  to  an  adequate  degree,  the 
meteorologically  important  part  of  the  forecast  can  be  permanently 
and  adversely  altered  by  the  model’s  response  to  the  initial  imbalance. 

Incompatibility  between  the  mass  and  momentum  fields  can  be 
generated  in  many  ways: 

(1)  errors  are  contained  in  the  observations, 

(2)  observations  have  contained  information  on  scales  of 
motion  not  represented, 

(3)  interpolations,  both  vertical  and  horizontal,  have 
reduced  the  balance,  and 

(4)  a  balanced  state  in  the  atmosphere  does  not  necessarily 
correspond  to  a  balanced  state  in  the  model  because  of 
the  numerous  approximations  contained  in  the  model. 

As  an  example  of  (4),  Nitta  and  Hovermale  (1967)  pointed  out  that 
boundary  conditions  and  finite  difference  methods  employed  in  the 
numerical  formulation  of  the  model  atmosphere  can  produce  gravity- 
inertia  waves. 

1.1.3  Balancing  of  the  mass  and  momentum  fields  in  PE  models 

The  purpose  of  balancing  is  to  reduce  the  generation  of 
nonmeteorological  gravity-inertia  waves  without  altering  the 
meteorologically  important  motions.  Unfortunately,  the  balanced  state 
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for  complicated  models  that  include  friction,  cumulus  convection, 
radiation,  etc. ,  is  usually  unknown. 

There  are  two  common  and  several  less  common  types  of  balancing 
methods  (Bengtsson,  1975;  Baer,  1977).  The  types  of  balancing  methods 
are: 

(1)  static  initialization  (SI), 

(2)  dynamic  initialization  (DI) , 

(3)  initialization  by  statistical  methods, 

(4)  initialization  utilizing  dynamical  integral  constraints,  and 

(5)  normal  mode  initialization. 

They  can  be  summarized  as  follows: 

(1)  SI.  In  a  static  initialization  (SI),  information 
applicable  primarily  at  one  time  is  used.  Time  derivatives  in  the 
equation  of  motion  are  usually  neglected  to  get  a  diagnostic  approxima¬ 
tion  to  the  balanced  state.  Examples  of  diagnostic  relationships  that 
have  been  used  are  the  geostrophic  wind  equation,  the  gradient  wind 
equation,  and  the  balance  equation. 

(2)  DI.  The  other  common  balancing  method  has  been  dynamic 
initialization  (DI).  In  DI,  the  actual  model  equations  are  used 
complete  with  time  derivatives.  A  DI  may  utilize  data  from  one  time 
only.  However,  some  DI  techniques  employ  data  from  two  or  more  times. 

In  this  type  of  DI,  sometimes  called  four-dimensional  data  assimilation, 
a  preforecast  integration,  usually  forward  in  time,  is  performed. 

(3)  Initialization  by  statistical  methods.  The  National 
Meteorological  Center  (NMC)  of  the  National  Weather  Service  (NWS)  has 


used  no  separate  initialization  for  their  operational  6-level  PE 
model.  The  forecast  is  started  directly  from  the  products  of  the 
objective  analysis.  The  model  is  initialized  with  analyzed  geopotentials, 
analyzed  nondivergent  winds,  and  12-hour  forecast  divergent  winds.  This 
procedure  was  adopted  because  the  short-range  forecast  errors  were 
smaller  than  when  the  balance  equation  was  used. 

(4)  Initialization  with  dynamic  integral  constraints.  To 
minimize  noise,  the  initial  analyses  can  be  adjusted  to  each  other  such 
that  they  satisfy  some  given  equation(s).  The  most  successful  approach 
has  been  that  of  Sasaki  (1958)  which  is  based  on  the  calculus  of 
variations.  First,  for  each  meteorological  variable,  the  difference 
between  adjusted  and  observed  quantities  at  a  point  is  expressed  as  a 
sum  of  squares.  Then  the  integral  over  the  volume  of  all  the  sums  of 
squares  is  minimized  such  that  a  given  constraint  is  satisfied.  This 
method's  major  drawback  has  been  its  complexity. 

(5)  Normal  mode  initialization.  First,  the  normal  inodes  of 
response  for  the  model  to  be  initialized  must  be  determined.  For 
example,  a  simple  baroclinic  model  may  have  one  fast  external  gravity- 
inertia  wave  mode  and  one  slow  internal  gravity-inertia  wave  mode. 

The  observations  are  then  resolved  into  a  series  expansion  of  the 
normal  modes  in  such  a  way  as  to  detect  the  fast  and  slow  modes. 

Then  the  coefficients  of  the  fast  modes  are  set  equal  to  zero.  With 
this  method,  the  initial  vcrticity  field  is  presumed  correct  and  the 
corresponding  geopotential  and  divergence  fields  that  would  eliminate 
high-frequency  noise  are  found.  This  initialization  technique  is 
extremely  complex,  time  consuming,  and  is  just  beginning  to  be  used. 
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1.2  Previous  research  on  the  initialization  of  PE  models 

1.2.1  Historical  review  of  initialization 

The  early  work  on  initialization  consisted  primarily  of  replacing 
subjective  techniques  for  the  analysis  of  observations  to  grid  points 
by  objective  analysis.  The  first  objective  initialization  technique 
was  due  to  Panofsky  (1949).  Panofsky  approximated  wind  and  pressure 
observations  by  third-degree  polynomials  requiring  10  coefficients. 
Boundaries  and  data-sparse  regions  were  difficulties  encountered  with 
this  method. 

The  successive  correction  type  of  objective  analysis  was  first 
introduced  by  Bergthdrsson  and  DSSs  (1955).  First-guess  values  of 
geopotential  were  obtained.  These  values  were  then  corrected  by 
geopotential  and  wind  observations  in  one  additional  scan  over  the  grid. 

Cressman  (1959)  also  employed  the  successive  correction  method  of 
objective  analysis.  He  made  several  passes  with  an  increasingly  small 
radius  of  influence.  Because  of  the  method's  simplicity  and  its 
applicability  to  most  atmospheric  variables,  even  in  data-sparse 
regions,  it  was  the  primary  analysis  tool  of  the  NWS  until  its 
replacement  in  1974. 

As  models  became  more  sophisticated,  research  on  balancing 
became  more  important.  One  of  the  earliest  studies  of  mass-momentum 
balancing  was  by  Hinklemann  (1951).  Using  a  linear  barotropic  model, 
he  showed  that  the  amplitude  of  the  undesired  gravity-inertia  waves 
was  reduced  by  replacing  the  observed  winds  with  geostrophic  winds. 
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Charney  (1955)  utilized  the  balance  equation  to  obtain  geopotentials 
from  the  wind  field.  Using  a  primitive  equation  model  and  artificial 
winds,  he  showed  that  the  generation  of  large  amplitude  gravity-inertia 
waves  was  greatly  suppressed  when  geostrophic  winds  were  replaced  by 
balanced  winds. 

1.2.2  Dynamic  initialization 

Dynamic  initialization  was  first  studied  by  Nitta  and  Hovermale 
(1967).  To  balance  the  mass  and  momentum  fields,  the  initial  conditions 
on  mass  and  momentum  were  forecast  one  time  step  forward  followed  by 
one  time  step  backward.  Then  the  mass  or  momentum  field  was  restored 
to  its  initial  value.  The  cycle  was  repeated  numerous  times.  The 
resulting  gravity-inertia  waves  were  damped  by  the  Matsuno  (1966)  time 
differencing  scheme.  It  was  shown  that  the  model  equations  could  be 
used  to  achieve  balanced  initial  conditions. 

A  dynamic  initialization  technique  was  also  studied  by  Hoke  and 
Anthes  (1976).  In  this  technique,  a  12-hour  pre-forecast  integration 
is  performed  in  which  the  model  variables  are  nudged  toward  the  observed 
value  of  the  variables  at  each  grid  point  and  time  step.  Subsequent 
12-hour  forecasts  were  better  than  those  forecasts  started  from  a  SI. 
However,  this  kind  of  DI  had  disadvantages.  First,  it  required 
additional  computer  time  (in  this  case  24  hours  of  simulated  time  for 
a  12-hour  forecast).  Second,  at  the  end  of  the  pre-forecast  integration, 
the  model  variables  have  not  necessarily  been  close  to  the  known  grid 


point  values. 


Temperton  (1976)  performed  experiments  in  dynamic  initialization 
with  a  5-level  hemispheric  model.  The  model  was  run  10  days  to  achieve 
a  balanced  state  which  was  considered  the  control  run.  Temperton 
found  that  external  gravity  waves  led  to  smaller  forecast  rain  amounts 
than  occurred  without  external  gravity  waves.  Additionally,  he  found 
that  dynamic  initialization  yielded  smaller  rain  amounts  than  the 
control.  Temperton  showed  that  external  gravity  waves  should  be 
rapidly  damped  by  DI  while  the  low  frequency  internal  modes  are 
relatively  unaffected. 

1.2.3  Normal  mode  analysis  and  initialization 

Phillips  (1960)  and  Blumen  (1975)  have  suggested  incorporating 
the  model's  normal  modes  of  response  into  the  initialization  process. 
Flattery  (1970)  employed  Hough  functions  in  his  analysis  scheme  in 
which  atmospheric  data  are  expanded  in  terms  of  the  Rossby  modes  of 
free  oscillation  of  the  shallow  fluid  equations.  This  implied  that 
no  gravity-inertia  waves  should  be  generated  in  an  exact  integration 
of  the  linear  equations. 

Dickinson  and  Williamson  (1972)  suggested  linear  normal  mode 
initialization  and  determined  the  normal  modes  of  a  2-level  model 
based  on  the  shallow  fluid  equations  expanded  in  terms  of  spherical 
harmonics.  Williamson  (1976)  then  applied  a  normal  mode  initialization 
to  a  shallow  water  grid-point  model.  He  found  reduced  gravity-inertia 
oscillations  in  nonlinear  integrations  but  the  gravity-inertia  waves 


were  not  eliminated. 
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Williamson  and  Dickinson  (1976)  determined  the  normal  response 
modes  (vertical  and  latitudinal)  for  the  linearized  National  Center 
for  Atmospheric  Research  (NCAR)  general  circulation  model  (GCM) . 

The  ensuing  integrations  were  found  to  contain  less  noise. 

Machenhauer  (1977)  and  Baer  (1977)  have  extended  normal  mode 
initialization  to  include  nonlinearity.  Machenhauer  first  determined 
the  free  normal  modes  of  a  spectral,  hemispheric,  shallow  fluid  model. 
He  then  determined  which  nonlinear  interactions  between  normal  modes 
lead  to  gravity-inertia  waves.  The  time  derivative  for  the  normal 
mode  coefficients  for  unwanted  modes  was  set  to  zero.  Machenhauer 
also  stated  that  normal  mode  initialization  produced  significant 
changes  in  the  forecasts  and  is  complicated  and  time  consuming. 

Another  disadvantage  of  this  type  of  initialization  has  been  that 
if  the  model  is  changed,  the  normal  response  modes  are  changed,  and 
the  model  normal  modes  must  be  reanalyzed. 

Baer  (1977)  applied  normal  mode  initialization  by  allowing 
nonlinearity  to  affect  the  initial  conditions  in  such  a  way  as  to 
eliminate  gravity- inertia  waves  in  the  initial  conditions  and  prevent 
most  of  them  from  occurring  during  the  time  integration.  Baer  assumed 
the  initial  vorticity  was  accurate.  He  then  adjusted  the  geopotentials 
and  divergence  so  they  were  compatible  with  the  vorticity  field.  In 
effect,  the  high-frequency  waves  (external  gravity  and  fastest 
internal  gravity  waves)  were  eliminated.  However,  the  slower  gravity- 
inertia  modes  have  speeds  comparable  to  Rossby  modes  and  were  not 


affected. 


15 


Baer  and  Tribbia  (1977)  extended  Baer's  (1977)  technique  to 
any  prediction  model  of  a  planetary  fluid  of  reasonably  small  Rossby 
number. 


1.2.4  Geostrophic  adjustment 

Geostrophic  adjustment  was  discussed  briefly  in  Section  1.1.2. 
However,  0kland  (1970)  reached  some  conclusions  which  are  important 
to  this  work  and  they  will  be  presented  here. 

0kland  found  a  linear  solution  to  a  simple  linear  baroclinic 
model  for  given  initial  conditions.  The  solution  consisted  of  two 
parts:  (1)  high-frequency  gravity-inertia  waves,  and  (2)  low-frequency 

gravity- inertia  waves  which  effect  a  balance  between  the  mass  and 
momentum  fields.  0kland  believed  that  the  solution  to  the  general 
analogous  initial-value  problem,  if  it  could  be  derived,  would  consist 
of  two  parts;  the  one  being  high-frequency  gravity- inertia  waves  and 
the  other  the  balanced  mass  and  wind  fields  which  exhibit  much  lower 
frequency. 

0kland  stated  that  gravity  waves  may  be  suppressed  by  damping  or 
propagation  away  from  the  area  of  interest.  A  damping  integration 
scheme  such  as  the  Euler-backward  scheme  could  be  used  but  the  slower 
internal  gravity  waves  move  at  about  the  same  speed  as  the  Rossby  modes. 
Hence  damping  the  slower  gravity-inertia  waves  would  also  damp  the 
waves  of  interest.  Dispersion  of  the  wave  energy  from  the  source  would 
be  effective  for  waves  with  high  group  velocities.  If  the  initial 


data  created  waves  of  small  group  velocity,  then  those  waves  would 


disperse  slowly  and  the  adjustment  would  be  slow,  a  situation  to  be 
avoided.  In  contrast  with  baroclinic  models,  barotropic  models  of 
deep  motions  contain  only  one  mode  and  it  is  fast  (an  external  gravity 
wave  with  a  speed  of  about  300  m  s  ^).  Therefore,  the  adjustment  in 
barotropic  models  is  relatively  rapid. 

0kland  conducted  two  experiments  with  a  2-level  baroclinic  model. 
The  model  was  allowed  to  run  for  24  hours.  Experiment  I  was  the  next 
24  hours  of  the  forecast.  In  Experiment  II,  the  divergence  was 
removed  from  the  initial  conditions  before  the  second  24-hour  forecast 
was  begun.  He  plotted  the  root-mean-square  (RMS)  values  of  the  local 
rate  of  change  of  surface  pressure  and  the  RMS  values  of  omega  (the 
vertical  velocity  in  pressure  coordinates)  at  about  500  mb.  The 
pressure  curves  were  more  indicative  of  external  gravity  wave  activity. 
Gravity  wave  activity  can  be  visualized  as  the  effect  of  all  gravity 
waves  simultaneously  acting  to  modify  the  balance  between  the  mass  and 
momentum  fields.  This  activity  was  greatest  for  the  nondivergent 
initial  conditions  (Experiment  II).  The  RMS  omega  graphs  were  more 
indicative  of  the  internal  gravity- inertia  wave  activity.  This 
activity  contained  less  noise  in  the  case  that  contained  the  divergence 
(Experiment  I). 

1.2.5  Balancing  on  sigma  versus  pressure  surfaces 

Most  PE  models  today  employ  pressure  or  normalized  pressure  (sigma) 
as  the  vertical  coordinate.  Reasons  include:  upper  air  data  are 
reported  on  standard  pressure  levels;  density  does  not  appear  explicitly 
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in  the  model  equations;  and,  the  boundary  conditions  at  the  ground  are 
easier  to  formulate.  The  mesoscale  model  developed  at  Penn  State  (see, 
e.g.,  Anthes  and  Warner,  1978,  hereafter  referred  to  as  AW)  uses  the 
sigma  coordinate  system.  In  the  initialization  of  this  model, 
geopotentials  are  calculated  on  pressure  surfaces.  Balanced  temperatures 
are  then  derived  and  the  winds  and  temperatures  are  interpolated  to  the 
model  sigma  surfaces. 

Since  the  nonlinear  balance  equation  has  been  used,  the  model 
initial  conditions  have  been  nondivergent  and  the  initial  vertical 
velocity  (omega)  at  each  grid  point  should  have  been  zero.  However, 
Warner  et  al.  (1978,  hereafter  referred  to  as  WAM)  reported  that  small 
omega  values  on  the  order  of  10  mb  d  ^  are  introduced  on  the  sigma 
surfaces  in  the  interpolation  from  pressure  to  sigma  surfaces.  In 
other  words,  the  initial  balance  has  been  altered.  These  omega  values 
on  sigma  surfaces  are  largest  where  terrain  slopes  are  greatest  and 
hence  have  not  necessarily  corresponded  to  meteorological  features. 

Since  the  model-generated  omega  values  were  much  larger,  the  small 
initial  omega  values  have  not  been  considered  to  be  a  problem. 

An  alternative  approach  has  been  to  interpolate  the  model 
variables  to  sigma  surfaces  prior  to  balancing.  Sundqvist  (1975) 
interpolated  observed  geopotentials  to  sigma  surfaces  and  then  used 
the  nonlinear  balance  equation  to  calculate  the  stream  function  on 
sigma  surfaces.  He  used  a  hemispheric  five-level  model  with  a  grid 
length  of  300  km  at  60°N.  His  scale  analysis  of  the  divergence  equation 
in  sigma  coordinates  indicated  that  the  nonlinear  balance  equation  was 
appropriate  for  his  model.  The  initialization  on  sigma  surfaces 
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reduced  initial  gravity- inertia  wave  oscillations  when  compared  with 
initialization  on  pressure  surfaces  and  interpolation  to  sigma 
surfaces.  However,  an  initialization  on  sigma  surfaces  with  orography 
included  produced  more  gravity- inertia  wave  oscillations  than  an 
initialization  with  no  terrain. 

1.2.6  Effect  of  divergent  initial  conditions 

Phillips  (1960)  pointed  out  that  the  removal  of  gravity-inertia 
waves  from  even  simple  models  is  not  just  a  matter  of  specifying 
initial  geostrophic  or  nondivergent  flow.  Small  values  of  divergence 
were  actually  needed  in  the  initial  conditions  to  suppress  these 
waves.  To  specify  the  initial  divergence  in  baroclinic  models, 

Phillips  proposed  various  forms  of  the  geostrophic  omega  equation.  The 
adiabatic  geostrophic  omega  equation  relates  vertical  motion  to  the 
advection  of  temperature  and  vorticity. 

Warner  (1972)  initialized  jets  in  barotropic  channel  and 
hemispheric  models.  He  used  the  following  initializations: 

(1)  geostrophic, 

(2)  mass  field  in  balance  with  the  fully  divergent  wind  field, 

(3)  mass  field  in  balance  with  the  nondivergent  wind  field, 

(4)  quasi-gradient,  and 

(5)  backward-forward  integration  about  the  initial  time  (a 
type  of  DI). 
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The  geostrophic  initialization  was  clearly  inferior  as  curvature  is 
neglected.  The  most  important  finding  relevant  here  was  that  the 
initializations  including  divergence  in  the  initial  conditions  ((2) 
and  (5))  produced  the  least  gravity-inertia  wave  noise. 


1.2.7  Previous  divergent  initialization  techniques  and  their  effect 
on  precipitation  predictions 


The  initialization  of  the  divergent  component  (hereafter  referred 
to  as  divergent  initialization)  has  been  attempted  by  several 
researchers.  Houghton  et  al.  (1971)  obtained  omega  values  from  a 
diagnostic  equation  similar  to  the  omega  equation.  Their  omega 
equation  did  not  contain  a  diabatic  term.  The  coarse-mesh  grid  used 
had  grid  points  at  five  degree  intervals  in  latitude  and  longitude.  The 
divergent  initial  conditions  had  a  small  effect  on  the  dynamic  variables 
at  first,  no  effect  after  12  hours,  and  the  same  level  of  noise  as  the 
nondivergent  initial  conditions.  The  initial  conditions  were  not 
saturated  where  vertical  motion  was  present  and  the  effect  of  divergent 
initialization  on  precipitation  was  not  evaluated.  They  stated  that 
it  may  be  that  poor  model  resolution  of  the  small  scales  inherent  in 
the  vertical  motion  field  contributed  significantly  to  the  lack  of 
forecast  improvement. 

LejenMs  (1977)  used  the  quasi-geostrophic  omega  equation  without 
a  diabatic  term  to  obtain  omega  fields  on  a  coarse-mesh  grid  (grid 
increment  300  km).  The  continuity  equation  was  then  used  along  with 
a  relaxation  procedure  to  convert  the  omega  values  to  velocity 


20 


potentials.  After  adding  the  divergent  component  to  the  initial 
conditions,  he  found  a  higher  precipitation  rate  than  the  nondivergent 
forecast  for  the  first  four  hours,  a  lower  precipitation  rate  from 
4  to  10  hours,  and  the  same  noise  level  as  nondivergent  initial 
conditions.  The  moisture  field  initially  was  not  saturated.  The 
effect  of  divergent  initialization  on  the  dynamic  variables  was  not 
reported.  Although  the  precipitation  forecast  was  more  realistic,  it 
was  significantly  less  than  the  observed  precipitation.  We  would 
expect  that  after  the  model  has  reached  a  balanced  state  internally, 
the  initial  divergent  component  would  no  longer  be  important.  Instead 
the  model  would  have  produced  a  divergent  component  compatible  with 
Che  model's  nondivergent  component. 

Dey  and  McPherson  (1977)  initialized  the  divergent  component  in 
the  NMC  coarse-mesh  hemispheric  PE  model.  The  initial  balanced  state 
was  derived  from  observations.  The  divergent  component  was  derived 
from  a  model  forecast  valid  at  the  same  time.  A  vertical  velocity 
equation  was  not  used.  Dey  and  McPherson  had  thought  that  divergent 
initialization  might  be  beneficial  if  applied  over  several  forecast 
cycles,  but  the  divergent  initialization  caused  only  small  changes  in 
both  the  global  analyses  and  forecasts.  The  only  differences  in  the 
precipitation  forecasts  occurred  in  areas  of  very  light  precipitation 
and  therefore  the  significance  was  difficult  to  assess.  They  concluded 
that  divergent  initialization  neither  degraded  nor  improved  the 
performance  of  the  NMC  global  system.  Additionally ,  Dey  (personal 
communication)  reported  that  a  forecast  divergent  component  that  was 
applicable  5  days  before  the  initialization  time  was  inadvertently 
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added  to  balanced  initial  conditions  in  one  experiment  and  no 
significant  change  in  the  forecast  resulted. 

Lubeck  et  al.  (1977)  have  performed  divergent  initialization 
experiments  with  a  global  spectral  model.  They  concluded  that  the 
divergent  initialization  had  a  relatively  small  effect  on  the  forecast 
dynamic  variables.  A  semi-implicit  integration  scheme  was  used  to 
smooth  the  high  frequency  oscillations.  Orography  was  not  included 
in  the  experiments.  Also,  the  model  was  dry  and  therefore  the  effect 
of  the  divergent  initialization  on  precipitation  could  not  be  evaluated. 

Smagorinsky  et  al.  (1967)  did  experiments  with  a  global  9-level 
PE  model.  They  used  a  form  of  the  omega  equation  with  no  diabatic 
term.  The  experiments  started  with  a  divergent  initialization  produced 
almost  the  same  precipitation  as  forecasts  started  from  nondivergent 
initial  conditions. 

Benwell  et  al.  (1971)  used  a  10-level  PE  model  for  divergent 
initialization  experiments.  They  reported  slightly  larger  precipitation 
amounts  with  nondivergent  initial  conditions. 

Divergent  initialization  has  also  been  used  in  hurricane  models. 

In  hurricanes,  the  divergent  part  of  the  wind  is  large  and  forecasts 
beginning  with  the  total  wind  have  been  more  successful  than  those 
started  from  the  nondivergent  component  only.  Miller  et  al.  (1972) 
obtained  the  stream  function  and  then  temperatures  were  derived.  The 
quasi-geostrophic  omega  equation  with  a  diabatic  term  was  solved  for 
omega  values  which  are  in  turn  used  in  the  continuity  equation  to  obtain 
velocity  potentials  (presumably  the  relaxation  boundary  condition  was 
one  of  constant  velocity  potential).  The  total  wind  was  obtained  by 
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combining  the  divergent  and  nondivergent  parts.  Temperatures  and  winds 
were  then  forecast.  From  that  forecast,  the  time  dependent  terms  in 
the  complete  divergence  and  omega  equations  could  be  evaluated.  The 
heights,  temperatures,  and  omega  values  were  recalculated.  The  process 
was  repeated  until  the  heights  and  omega  fields  were  relatively  stable. 
This  process  was  similar  to  DI.  The  areal  extent  of  their  precipitation 
forecasts  agreed  with  the  observed  rainfall.  However,  the  forecast 
precipitation  amounts  were  too  low.  The  model  atmosphere  was  not 
saturated  initially  and  the  forecast  model  did  not  contain  a  convective 
precipitation  parameterization.  The  adequacy  of  the  initial  precipitation 
rate  was  not  evaluated. 

Mathur  (1974)  used  an  omega  equation  similar  to  the  quasi-geostrophic 
omega  equation  but  it  did  not  contain  a  diabatic  term.  The  omega  values 
were  then  used  in  the  continuity  equation  to  obtain  velocity  potentials. 
Again,  presumably  a  constant  value  was  used  as  the  boundary  condition 
on  the  velocity  potential.  The  meshed  grid  increments  were  37  km  on  the 
interior  and  74  km  on  the  exterior.  The  hurricane  studied  occurred 
entirely  over  the  ocean  and  observed  precipitation  fields  were  not 
available. 

In  summary,  there  has  been  a  lot  of  work  on  the  various  aspects 
of  the  initialization  problem.  The  work  to  date  on  divergent 
initialization  has  been  primarily  on  the  large  or  synoptic  scale.  On 
that  scale,  the  divergent  initialization  has  had  little  if  any  effect 
on  the  forecast.  The  divergent  initializations  on  the  mesoscale  have 
been  for  hurricane  models.  The  resulting  forecasts  were  evaluated  by 
comparing  forecast  versus  observed  track  and  intensification.  The 
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adequacy  and  suitability  of  the  divergent  initializations  as  compared 
with  nondivergent  initializations  were  not  evaluated.  The  effect  of 
divergent  initialization  on  precipitation  forecasts  was  not  closely 
scrutinized.  Finally,  diabatic  heating  information  and  asynoptic 
data  were  not  utilized  in  the  initialization  schemes. 

1.3  Synoptic  scale  versus  mesoscale  divergence  values 

Before  any  further  discussion  of  the  differences  between  the 
synoptic  and  mesoscales,  we  will  present  a  more  quantitative  definition 
of  the  scales.  In  Table  1,  we  adopt  Orlanski's  (1975)  scale  definitions 
but  we  have  subdivided  the  meso-B  scale  into  upper  and  lower  portions. 

In  subsequent  sections,  synoptic  scale  will  refer  to  meso-a  and  larger 
scales  and  the  term  mesoscale  will  refer  to  the  upper  meso-B  scale 
unless  otherwise  noted. 

Table  1.  Scale  Definitions. 

Scale  designation  Scale  range  (km) 


meso-y 

2-20 

meso-6 

20-200 

lower 

meso-B 

20-60 

upper 

meso-8 

60-200 

tneso-a 

200-2000 

macro-8 


2000-10000 


As  discussed  in  Section  1.2,  divergent  initialization  on  the 
synoptic  scale  had  little  if  any  influence  on  the  precipitation 
forecasts.  We  will  now  attempt  to  gain  some  insight  into  why  divergent 
initialization  has  not  been  successful  on  the  synoptic  scale  but  might 
be  successful  on  the  mesoscale  and  smaller  scales.  Let  V  be  a 
characteristic  velocity  and  Ax  a  characteristic  horizontal  model  grid 
increment.  For  the  scales  in  Table  1,  10  m  s  ^  is  a  characteristic 
velocity.  For  characteristic  grid  increments,  let  us  compare  a 
synoptic  model  grid  increment  Ax^  =  300  km,  an  upper  meso-8  scale  grid 
increment  Ax^  =  100  km,  and  a  lower  meso-g  scale  grid  increment 
Ax^  =  30  km.  We  can  now  compare  typical  divergence  values  for  the 
three  different  grid  increments. 

We  know  from  the  divergence  theorem  that  the  net  divergence  over 
the  global  500-mb  surface  is  zero.  We  would  expect  the  same  result  to 
be  approximately  correct  for  a  hemisphere.  In  fact,  as  we  will  discuss 
in  more  detail  in  Chapter  2,  the  net  divergence  on  a  constant  pressure 
surface  over  an  area  of  about  3500  km  on  a  side  averages  almost  to 
zero.  However,  at  grid  increments  on  the  order  of  Ax^,  there  exist 
phenomena  such  as  thunderstorms  that  possess  large  divergences.  These 
divergence  patterns  could  easily  dominate  over  a  small  grid  and  produce 
large  net  divergence  at  a  given  level. 

Since  V  is  the  same  for  each  length  scale,  and  since  divergence 
is  of  order  V/ Ax,  we  expect  an  order  of  magnitude  larger  divergence 
values  as  we  progress  from  Ax^  to  Ax^.  The  grid  increments  Ax^, 
and  Ax^  each  differ  by  a  factor  of  three  but  the  area  contained  within 
each  grid  square  is  almost  an  order  of  magnitude  different.  That  is. 
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C  O  /  9  *)  0 

(Ax^)  'v  10  km  ,  (Ax^)  ^  10  km  ,  and  (Ax^)  ^10  km  .  Therefore, 

2 

we  can  visualize  a  much  larger  net  divergence  over  an  area  (Ax^)  than 
over  an  area  (Ax^)^. 

Thus,  the  divergence  at  a  grid  point  can  be  larger  in  mesoscale 
models  than  in  synoptic  scale  models.  Because  divergence  is  the  forcing 
function  in  the  Poisson  equation  for  velocity  potential  (and  hence  the 
divergent  component),  we  expect  a  mesoscale  divergent  component  of 
greater  magnitude  than  on  the  synoptic  scale.  In  other  words,  the 
divergent  component  may  well  be  a  significant  ingredient  in  mesoscale 
model  initial  conditions  while  previous  research  has  indicated  that  it 
is  not  significant  in  that  respect  on  the  synoptic  scale.  We  will 
accomplish  a  scale  analysis  in  Chapter  2  to  see  how  significant  the 
divergence  can  be  expected  to  be  on  the  mesoscale  but  we  have  demonstrated 
here  in  qualitative  terms  that,  if  divergent  initialization  is  to 
significantly  affect  a  precipitation  forecast,  it  would  most  likely 
occur  on  the  mesoscale  and  smaller  scales. 

1.4  Research  objectives 

To  date,  precipitation  forecasting  with  NWP  models  has  not  been 
very  accurate,  especially  in  cases  of  significant  precipitation. 

Accurate  and  timely  precipitation  forecasts  on  the  mesoscale  would  be 
of  great  economic  and  other  value  in  times  of  flash  floods,  squall 
lines,  heavy  precipitation  over  large  areas  such  as  hurricanes,  etc. 
Therefore,  a  divergent  initialization  procedure  which  improves 
mesoscale  NWP  model  precipitation  forecasts  would  be  an  important 


contribution. 
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A  general  divergent  initialization  procedure  developed  in  this 
thesis  can  be  described  as  follows:  a  vertical  velocity  equation 
diagnoses  vertical  velocities  from  observed  meteorological  fields. 
Velocity  potentials  are  then  derived  from  the  vertical  velocities  with 
the  continuity  equation  and  knowledge  of  appropriate  boundary  conditions 
on  velocity  potential.  Next,  the  divergent  wind  components  are 
obtained  from  the  velocity  potential.  Geopotentials  are  calculated 
on  sigma  surfaces  using  a  balance  equation  with  contributions  from 
both  the  nondivergent  and  divergent  wind  components.  Finally,  balanced 
temperatures  are  derived  via  the  hydrostatic  equation.  Note  that  the 
vertical  velocity  and  divergence  fields  are  diagnosed.  We  cannot 
accurately  measure  divergence  directly  because  the  errors  in  the  wind 
observations  are  on  the  order  of  magnitude  of  the  divergent  wind 
component.  Therefore,  we  derive  a  vertical  velocity  equation  that 
defines  what  vertical  velocity  field  would  exist  for  a  given  set  of 
dynamic  and  thermodynamic  variables. 

Questions  to  be  addressed  in  this  thesis  include: 

General : 

(1)  What  are  the  differences  between  divergent  initialization 
on  the  synoptic  and  mesoscales? 

(2)  Does  the  inclusion  of  a  divergent  component  on  the 
mesoscale  improve  the  ensuing  forecast,  especially 
precipitation  rate  and  amount?  Analogously,  does  the 
model  remember  the  divergence  as  introduced  by  the 


initialization? 
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(3)  Do  modeled  precipitation  systems  recover  their  "correct" 
intensity  later  in  the  forecast,  after  a  nondivergent 
initialization,  or  are  the  precipitation  rates  throughout 
the  forecast  adversely  affected  by  the  incorrect  vertical 
velocities  and  latent  heat  release  produced  early  in  the 
forecast  by  the  lack  of  initial  divergence? 

(4)  What  form  of  the  omega  and  balance  equations  should  be 
used  on  the  mesoscale? 

Specific : 

(5)  Can  the  vertical  velocities  be  diagnosed  to  sufficient 
accuracy  to  be  useful  in  the  divergent  initialization 
procedure? 

(6)  What  boundary  conditions  should  be  used  on  stream  function, 
omega,  and  geqpotential? 

(7)  The  continuity  equation  is  used  to  derive  velocity 
potential  from  vertical  velocity.  How  should  the 
boundary  conditions  on  velocity  potential  be  specified 
to  achieve  maximum  accuracy  on  the  domain  interior? 

The  forecast  model  that  will  be  used  for  these  experiments  is 
the  mesoscale  NWP  model  developed  at  The  Pennsylvania  State  University 
(see  WAM,  AW).  In  Chapter  2  a  scale  analysis  of  the  relevant  equations 
as  well  as  a  thorough  analysis  of  the  appropriate  boundary  conditions 
for  the  required  dependent  variables  will  be  provided.  In  Chapter  3, 
a  brief  description  of  the  forecast  model  will  be  given.  The  data  at 
the  Initial  and  final  synoptic  times,  the  observed  precipitation,  and 
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a  discussion  of  the  synoptic  situation  chosen  for  study  will  be 
presented  in  Chapter  4.  In  Chapter  5,  we  will  develop  the  finite- 
difference  form  of  the  omega  equation  which  will  be  used  to  determine 
vertical  velocities.  The  omega  equation  will  then  be  applied  to  the 
data  set.  In  Chapter  6,  the  finite-difference  form  of  the  balance 
equation  will  be  derived  from  the  model  equations.  In  Chapter  7, 
the  finite-difference  balance  equation  will  be  applied  to  the  chosen 
data  set.  Chapter  8  will  contain  the  summary  and  conclusions. 
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2.0  DEVELOPMENT  OF  THE  DIVERGENT  INITIALIZATION  PROCEDURE 

In  this  chapter,  we  will  first  state  which  forces  are  considered 
in  the  forecast  model  developed  at  The  Pennsylvania  State  University 
(PSU).  Then  we  will  perform  a  scale  analysis  of  the  complete 
vertical  velocity  and  divergence  equations.  These  scale  analyses 
will  specify  which  terms  must  be  retained  in  the  respective  equations 
for  use  in  the  specific  divergent  initialization  procedure.  Finally, 
we  will  examine  which  boundary  conditions  should  be  employed  for  the 
various  second-order  elliptic  partial  differential  equations  (PDEs) 
which  arise  in  this  scheme.  We  shall  see  that  the  proper  specification 
of  boundary  conditions  is  important  to  this  initialization  procedure. 
The  general  framework  for  the  divergent  initialization  procedure  will 
be  outlined  whereas  the  specific  technique  for  the  initialization  of 
the  PSU  model  will  be  addressed  in  later  chapters. 

Before  we  can  perform  the  required  scale  analyses,  we  need  to 
know  what  forces  are  considered  by  the  PSU  model.  That  is,  we  need 
to  know  what  form  of  the  equation  of  motion  is  used  so  that  we  can 
derive  consistent  vertical  velocity  and  divergence  equations.  The  PSU 
model  considers  the  pressure  gradient  force,  the  Coriolis  force, 
gravity,  and  the  vertical  frictional  force  in  the  planetary  boundary 
layer  (PBL)  only.  Additionally,  horizontal  frictional  forces  are 
used  by  the  model  at  every  level,  but  for  numerical  rather  than 


physical  reasons. 
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2.1  Scale  analysis  of  the  vertical  velocity  equation 

Scale  analysis  is  a  technique  which  permits  the  determination 
of  what  terms  in  a  given  equation  are  important,  based  on  characteristic 
values  of  physical  variables  and  characteristic  scales  in  space  and 
time  of  the  class  of  phenomena  of  interest.  Here,  we  are  interested 
in  the  upper  meso-0  and  larger  scales  (e.g.,  baroclinic  waves). 

The  scale  analysis  presented  here,  through  (2.19),  closely 
follows  that  of  Williams  (personal  communication).  We  define 

L  =  characteristic  horizontal  scale  (roughly  a  quarter 
wavelength  of  the  disturbances  considered) 

T  *  characteristic  time  scale 

V  =  characteristic  horizontal  velocity 

.  The  approximate  magnitudes  of  derivatives  can  be  estimated  as  follows: 


3u  3v  V 
3x  3y  L 


(2.1a) 


3u 

at 


v 


(2.1b) 


The  time  scale  is  generally  given  by 


T 


(2.2) 


Here  we  will  use  V  'v  10  m  s  ^  for  the  synoptic  and  mesoscales,  L  %  1000  km 
for  synoptic-scale  motions  and  L  ^  100  km  for  mesoscale  features. 


The  Helmholtz  theorem  states  that  a  velocity  vector  can  be 
separated  into  nondivergent  and  divergent  (irrotational)  parts: 


V  =  V,  +  V 
-i)*  "X 


(2.3) 


where 


V  =  k  x  74i  (2. 3a) 


and 


V  =  VX  .  (2.3b) 

In  (2.3a)  and  (2.3b),  ip  Is  the  stream  function  and  x  is  the  velocity 
potential.  Since  we  know  that  the  magnitude  of  the  nondivergent  wind  is 
on  the  order  of  the  observed  wind  for  the  scales  (upper  meso-g  and 
larger)  considered,  we  use 


V,  ^  V 
i> 


(2.4) 


for  the  scale  of  the  nondivergent  wind.  On  the  other  hand,  the 
divergent  wind  is  generally  an  order  of  magnitude  smaller.  Therefore, 
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where  R^  ^  0.1,  but  for  the  upper  meso-8  and  larger  scales  of  motion 
which  we  consider  here,  could  be  as  small  as  zero  or  as  large  as 
0.3. 

In  addition  to  those  scaling  parameters,  we  will  use 


V  %  ~  ,  (2.6a) 

Ro  =  ?l  or  s‘hi  ’  (2-6b) 


D  ^  Rx  £  ,  (2.6c) 

C  -v  l  ,  (2 . 6d) 


and 


6  =  vf  =  -^  =  cos  $  4^  ^  ~  (2. 6e) 

dy  dy  a 

for  mid-latitudes.  If  the  wind  were  geostrophic,  then 


V<f>  -v  fV  .  (2. 6f ) 


Now  we  separate  the  geopotential  into  two  parts. 

4>  *  $( z)  +  <J>'  (2.6g) 
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In  (2.6g),  <J>'  is  the  perturbation  geopotential  and  <p(z)  is  taken  from 
the  standard  atmosphere.  With  (2.6a),  (2.6f),  and  (2.6g),  we  obtain 
(2. 6h) . 


ll 

L 


^  fV  or  <p'  ^  fVL 


(2. 6h) 


Even  though  the  wind  is  not  geostrophic  on  the  upper  meso-B  scale, 
(2.6h)  is  adequate  for  scaling  purposes. 

2.1.1  Pertinent  equations  in  Z  coordinates 


We  utilize  the  vertical  coordinate  Z  which  was  first  used  by 
Phillips  (1963): 


Z  =  -  ln(^-)  .  (2.7) 

o 

Z  is  related  to  actual  height,  h,  and  geopotential  tfi  by  the  equation  of 
state  and  the  hydrostatic  equation: 


RT  »  pa 


3Z 


From  (2.8)  we  see  that 


(2.8) 


i-  -  1- 

3p  "  3Z 


P 


(2.9) 


One  advantage  of  Z  is  that  it  approximately  equals  the  actual  height 
divided  by  the  scale  height.  That  is. 


Z 


(2.10) 


where  H  ■  RT/g  is  the  scale  height.  For  the  atmosphere,  H  <v  8  km. 
Therefore,  for  the  troposphere. 


Z  'v  1 


(2.10a) 


and 


3_ 

az 


1 


(2.10b) 


From  (2.7)  we  see 


(2.11) 


The  horizontal  vector  equation  of  motion  in  Z  coordinates  which 
contains  the  forces  considered  by  the  PSU  model  is 


av  .  av 

rf-  +  V-VV  +  Z  ry  +  V*  +  fk  x  V  -  Fr  *  0  .  (2.12) 

dt  ■*  ■*'  O  Lt  *  " 


The  continuity  equation  in  p  coordinates  is 
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Using  (2.11)  and  (2.9)  to  modify  (2.13),  we  obtain  the  continuity 
equation  in  Z  coordinates: 

3Z 

V-V  +  -  z  =  0  .  (2.14) 

With  the  hydrostatic  equation,  the  first  law  of  thermodynamics  can  be 
written 


If  f  +  i  - 2  If#  -  "♦>  ’ 


(2.15) 


where  <  =*  R/c  and  Q  is  the  heat  added  per  unit  time  and  mass. 
P 


2.1.2  Derivation  of  the  vertical  velocity  equation 

To  obtain  a  vertical  velocity  equation,  we  must  first  have  a 
vorticity  equation.  Operating  on  (2.12)  with  7  x  yields 


—■  +  V  -7?  +  V  *7?  +  V ,-Vf 

3t  ~ip  -X  ~'l» 


+  V  »7f  +  Z  rr-  +  fD  +  £D 
'X  8Z 


(2.16) 


ay  a  <  3v 

+  k*7Z  x  +  k-7Z  x  — -  7  x  Fr  =  0 


With  (2.6g),  (2.15)  can  be  rewritten  as 


36 


—  +  v  *v  +  v  *v  +  zr(z) 

at  az  -n»  az  -x  az  ^  ; 


+  z  fz(fz~  +  K*’)  ‘  kQ  =  0 


(2.17) 


where 


r(z)  =  az(H+  K*> 


R(f  +  *T) 


(2.18) 


JU (£_  +  I  il) 

T  Vc  H  3Z; 
T  P 


T  Is  a  static  stability  because  it  is  proportional  to  the  difference 
between  the  dry  adiabatic  lapse  rate  and  the  standard  atmospheric 
lapse  rate. 

We  derive  the  vertical  velocity  equation  by  operating  on  (2.16) 
3  2 

with  f  — ,  subtracting  V  of  (2.17),  and  reordering  terms: 
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(2.19) 
(cont . ) 


In  the  derivation  of  vertical  velocity  equations,  it  is  usually 
assumed  that  the  local  rate  of  change  of  the  actual  vorticity  equals 
the  time  rate  of  change  of  the  geostrophic  vorticity  (hereafter 
called  the  geostrophic  vorticity  assumption) ,  where  the  geostrophic 
vorticity  is 


Eg  •  7 


(2.20) 


We  make  this  assumption  here,  and  will  return  to  it  in  Section  2.1.5. 


2.1.3  Scale  analysis  of  the  vertical  velocity  equation 

From  the  continuity  equation  (2.14),  we  obtain 

Z  -v  RL  ^  (2.21) 

V 

since  D  m  R^  —  . 

Table  2  contains  a  scale  analysis  of  (2.19)  which  includes  a 

scaling  factor  for  each  forcing  function.  First,  consider  the  synoptic 

scale,  where  L  %  1000  km.  For  mid-latitudes,  that  implies  Ro  'v  0.1. 

The  value  of  is  taken  to  be  'v  0.001  (Anthes,  1978)  and  H  is  %  10  km. 

For  this  length  scale,  L/a  'v  0.1.  We  choose  a  precipitation  rate  of 

0.1  cm  h  1  as  typical  of  a  synoptic-scale  grid  square  (R2  =  0.1  in  T8  in 

Table  2).  Examination  of  the  scale  factors  of  (2.19)  reveals  that 

-14 

the  forcing  functions  T4,  T6,  T8,  and  Til  are  of  order  10  while 
terms  T5,  T7,  T9,  T10,  T12,  T13,  and  T14  are  an  order  of  magnitude 
smaller.  The  largest  forcing  functions  represent  differential 
vorticity  advection  by  the  nondivergent  wind,  the  Laplacian  of 
temperature  advection  by  the  nondivergent  wind,  the  diabatic  term,  and 
the  beta  or  Coriolis  term.  These  are  important  physical  processes  on 
the  large  scale.  For  divergent  initialization  on  the  large  scale,  the 
four  largest  forcing  functions  should  be  used  in  the  vertical  velocity 
equation.  It  should  be  noted  that  three  of  the  four  largest  terms 
(T4,  T6,  and  Til)  are  the  forcing  functions  of  the  quasi-geostrophic 
vertical  velocity  equation.  This  is  a  well-known  result  (e.g.,  see 
Haltiner,  1971). 
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Table  2.  Scale  analysis  of  the  forcing  functions 
in  the  vertical  velocity  equation. 


Term 

designator 


Scale 

factor 


T4 


T5 


TULVV1V,LV: 

Ro  L  L  L  Ro  l3 
^1  V3 

T5"R1 


T6 


T7 


1  IV 

T6  %  Vf  VL  -u  ■=— 

L2  R°  L3 

R1  V3 

T7"R1  T6"r^ 


T8 


Let  R.  be  a  precipitation  rate  of  1  cm  h  .  Then, 

^  2 
for  each  cm  , 


1  g  h 


a-  600  cal  h 


^  2400  joule  h 
"V  1  joule  s  3 


-1 


Therefore,  Q  ^  R£  per  unit  mass  in  mks  units  and 


T9 


kR_  R-R- 
T8  'V  — y  ^  — ~7~ 

L  L 

T9*^R 

Ro  L  1  L  L  Ro  l3 
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Table  2  (Continued) 


Term 

designator 


Scale 

factor 


T10 


Til 

T12 

T13 

T14 

T15 

T16 


The  model  friction  term  in  the  bulk-PBL  parameterization 
(see  AW)  can  be  scaled 


Z£„lL.Iil„iL-l3_p  cv2 

P*  P*  0  32  P*  a  32  s  D 


%  JL.  I  c  v2  V. 

RT  a  D  H 


c  ^ 

1  V  1  1  2  1  D  V 

Therefore,  T!0  ^  r  ^  CDV  /  3 


2  3 

IV  1  f  V  1  L  V 

Ro  L  Ro  a  L  R2  a  ^3 

R1  L  V3 

T12  %  R.  Til  ^  ~  - 

1  Ro  3  L3 

D  'l 

1  V  V  V  TV 

T13^HRir^ 


T14  'v  T13  -v 


hyl 

Ro  l3 


R1  V3 

T15  'v  R.  T14  -v  -r- 
1  Ro  L3 


1  V  V  V  IV' 
T16  -v  —  —  —  —  %  —  — 
Ro  L  L  L  Ro 


VI  ,TrT  1  V  V  IV 

T17  ^  r  ~o  fVL  ^  n - T71-; - T 

L  T  2  Ro  T  2  L  Ro  ,3 
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2.1.4  Determination  of  the  mesoscale  vertical  velocity  equation 

On  the  mesoscale,  L  *\>  1000  km.  For  mid-latitudes,  that  implies 
Ro  ^  1  and  L/a  'v  0.01.  We  take  a  precipitation  rate  of  1  cm  h  ^  as 
being  typical  of  a  mesoscale  grid  square  (R^  =  1  in  Table  2) .  Using 
these  values,  the  relative  order  of  magnitude  of  the  forcing  functions 
in  (2.19)  except  for  T16  and  T17  is  indicated  in  Roman  numerals 
underneath  the  terms  in  (2.19).  The  largest  forcing  function  is  the 
diabatic  term,  T8 ,  with  a  value  of  order  10  Therefore,  in  areas 

where  precipitation  is  occurring,  the  diabatic  term  dominates  the  other 
forcing  functions  on  the  mesoscale.  Terms  with  a  la  (T4  and  T6)  are 
an  order  of  magnitude  smaller  in  precipitation  areas.  These  terms 
represent  differential  vorticity  advection  and  the  Laplacian  of 
temperature  advection.  Note  that  these  terms  will  be  the  most  important 
terms  in  areas  where  no  precipitation  is  occurring.  Terms ^with  a  II 
(T5,  T7,  T9,  T10,  T13,  and  T14)  are  of  the  same  order  but  are  an 
order  of  magnitude  smaller  than  the  terms  with  a  la.  Forcing 
functions  with  a  III  (Til,  T12,  and  T15)  are  of  the  same  order  but  are 
an  order  of  magnitude  smaller  than  the  terms  labeled  II. 

The  goal  of  this  research  is  to  develop  and  test  a  divergent 
initialization  technique  on  the  mesoscale.  Since  the  divergent 
component  is  an  order  of  magnitude  smaller  than  the  nondivergent 
component  on  the  mesoscale,  and  since  we  desire  to  diagnose  the 
vertical  velocities  as  accurately  as  possible,  we  will  retain  terms  of 
first  and  second  order  and  discard  only  those  terms  at  least  two 
orders  of  magnitude  smaller  than  the  most  significant  terms.  That  is, 
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we  will  keep  terms  labeled  I,  la,  and  II  while  discarding  terms 
labeled  III.  Therefore,  the  final  mesoscale  vertical  velocity 
equation  consists  of  terms  Tl,  T2,  T3a,  T4,  T5,  T6,  T7,  T8,  T9,  T10, 
T13,  and  T14  of  (2.20)  and  will  be  referred  to  as  (2.22). 

The  omega  equation  corresponding  to  (2.22)  is 
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where  a  is  the  static  stability, 
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Beneath  each  term  is  a  Roman  numeral  indicating  the  term's  relative 
importance  on  the  mesoscale.  Eq.  (2.23)  with  appropriate  boundary 
conditions  will  be  used  in  Chapter  5  to  diagnose  omega  fields  from 
real  data. 

It  is  interesting  to  note  the  differences  between  the  synoptic 
and  mesoscales.  The  beta  or  Coriolis  term  is  a  first-order  term  on 
the  synoptic  scale  but  is  negligible  on  the  mesoscale.  That  is,  as 
we  go  to  smaller  scales,  the  Coriolis  force  is  less  important  and 
the  wind  is  less  geostrophic.  The  other  significant  difference  between 
the  scales  is  that  although  the  diabatic  term  is  important  on  the 
synoptic  scale,  it  is  relatively  more  important  on  the  mesoscale. 

While  the  diabatic  term  should  be  included  on  the  synoptic  scale,  it 
is  even  more  important  to  include  it  on  the  mesoscale.  In  other  words, 
local  forcing  by  latent  heating  is  more  important  on  smaller  scales. 

VJe  expect  that  the  vertical  velocities  to  be  diagnosed  on  the  mesoscale 
in  Chapter  5  will  reflect  the  diabatic  effect  to  a  larger  degree  than 
the  other  forcing  functions. 

For  divergent  initialization  to  succeed,  the  Initial  divergence 
must  be  remembered  by  the  model.  That  is,  it  must  be  supported  by 
the  model.  Otherwise,  the  initial  divergence  will  be  dissipated  by 
internal  gravity  waves.  The  initial  divergence  must  be  balanced  by 
latent  heating,  especially  on  the  mesoscale.  If  the  divergence  is  not 
maintained  by  latent  heating,  it  will  not  be  remembered.  Latent 
heating  should  occur  immediately  (in  the  first  time  steps  of  the 
forecast)  in  areas  of  upward  motion.  Therefore,  those  areas  should 
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initially  be  saturated  for  stable  layers  and  the  convective  parameteriza¬ 
tion  scheme  should  produce  convective  heating  in  unstable  moist  layers. 

As  stated  in  Chapter  1,  divergent  initialization  on  the  synoptic 
scale  has  had  little  effect.  There  was  little  effect  on  the  dynamic 
variables  and  a  slight  if  any  effect  on  the  initial  precipitation 
rates.  Where  the  omega  equation  was  used  in  the  determination  of  the 
divergent  component,  the  omega  equation  did  not  contain  a  diabatic 
term.  Although  the  initial  moisture  fields  were  not  given,  we  speculate 
they  were  not  initially  saturated  in  areas  of  upward  motion.  Hence  the 
initial  divergence  would  not  be  sustained  by  the  release  of  latent 
heat.  Also,  some  of  the  terms  in  the  vertical  velocity  equation  that 
are  most  important  on  the  synoptic  scale  such  as  vorticity  advection 
are  more  susceptible  to  initialization-related  noise.  It  is  therefore 
possible  that  the  initial  divergence  was  dissipated  before  it  could 
effect  the  initial  precipitation  rate.  We  will  later  show  that  on  the 
mesoscale,  the  initial  divergence  is  retained  by  the  forecast  model. 

It  is  possible  that  the  divergent  initialization  procedure  presented 
here  could  produce  better  results  than  those  previously  obtained  if 
applied  on  the  synoptic  scale. 

2.1.5  The  geostrophic  vorticity  assumption 

In  the  derivation  of  the  vertical  velocity  equation  in  Section 
2.1.2,  we  assumed  that  the  local  temporal  rates  of  change  of  the 
actual  vorticity  and  the  geostrophic  vorticity  were  equal.  Without 
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this  assumption,  the  vertical  velocity  equation  would  be  predictive 
rather  than  diagnostic. 

It  is  well  known  that  the  geostrophic  vorticity  assumption  is  a 
good  approximation  on  the  synoptic  scale  where  quasi-geostrophic 
theory  applies.  In  fact,  this  assumption  has  often  been  used  to 
obtain  the  vorticity  field  from  the  observed  height  field  rather  than 
from  wind  observations.  Scale  analysis  was  applied  to  terms  T16  and 
T17  of  (2.19)  as  shown  in  Table  2.  We  see  immediately  that  the  terms 
are  almost  equal  in  magnitude  and  opposite  in  sign.  It  first  appears 
that  the  geostrophic  assumption  will  be  valid  at  all  scales.  We  know, 
however,  that  where  the  effects  of  curvature  become  important  or 
where  the  divergent  component  is  large  such  as  in  gravity  waves,  the 
geostrophic  vorticity  assumption  is  not  at  ail  accurate.  For  gravity 
waves,  the  scale  analysis  breaks  down  because  the  advective  time 
scale  is  no  longer  appropriate  for  terms  T16  and  T17  of  (2.19). 
Therefore,  if  this  divergent  initialization  scheme  were  applied  in 
the  future  to  smaller  scales,  terms  T16  and  T17  may  have  to  be  included 
in  the  vertical  velocity  equation  used. 

Finally,  on  length  scales  where  the  geostrophic  vorticity 
assumption  breaks  down,  the  terms  labeled  with  a  III  in  (2.19)  may  not 
be  negligible  in  comparison  with  the  I  and  II  terms.  In  other  words, 
the  specific  vertical  velocity  equation  used  for  the  divergent 
initialization  of  the  PSU  model  depends  on  the  smallest  length-scale 
feature  that  will  be  permitted.  The  smallest  permissible  length  scale 
depends  on  the  grid  increment,  the  data  availability,  the  method  of 
analysis,  and  the  degree  of  smoothing. 
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2.2  Scale  analysis  of  the  divergence  equation 

This  scale  analysis  will  be  conducted  analogously  to  the  one  in 
Section  2.1.  The  same  scaling  relationships  will  be  used. 

2.2.1  Derivation  of  the  divergence  equation 

Recall  that  (2.12)  is  the  equation  of  motion  in  Z  coordinates. 
We  obtain  the  divergence  equation  by  operating  on  (2.12)  with  7*  to 
obtain 


3V  .  9V 

v"at  +  v* t  (V- V)V]  +  v-(z  pp 


+  vV  +  V-(fk  x  y)  -  V*Fr  *  0 


Using  (2.3),  (2.24)  becomes 
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(2.25) 
(cont . ) 

-  k  x  Vf*V  -  V*  Fr  =  0 
-X 

T5c  T6 

III  III 

which  is  the  complete  divergence  equation  in  Z  coordinates. 

As  is  customarily  done,  we  will  neglect  the  time-dependent  term 
Tl.  This  leaves  a  diagnostic  relationship  which  can  be  termed  the 
complete  balance  equation.  In  Section  2.2.4,  we  will  analyze  the 
effect  of  neglecting  this  term. 

2.2.2  Scale  analysis  of  the  balance  equation 

We  will  examine  (2.25),  term  by  term  (see  Table  3).  For  the 
synoptic  scale,  where  L  'v  1000  km  and  Ro  'v  0.1,  terms  T2a,  T4,  T5a, 
and  T5b  are  at  least  an  order  of  magnitude  larger  than  the  other  terms. 
These  terms  are  the  advection  of  momentum,  the  Laplacian  of  geopotential, 
the  Coriolis  parameter  times  the  relative  vorticity,  and  the  beta  term, 
respectively.  The  equation  containing  only  these  terms  is  commonly 
referred  to  as  the  nonlinear  balance  equation.  A  form  of  that 
equation  has  been  used  in  the  nondivergent  initialization  of  the  PSU 
model  (WAM;  Keyser,  1978;  Anthes,  1978).  It  has  been  previously 
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Table  3.  Scale  analysis  of  the  divergence  equation. 


Term 

designator 


Scale 

factor 


T1 


V  V  V 

T1  %  I  Ri  I  ^  Ri  T 
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T2a 


11  V 
T2a  *  -  V  -  V  ^  -j 

I. 
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T2b  ^  R  T2a  %  R.  -j 
1  1  L 
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T2c  ^  R.  T2a  ^  R.  —z 
1  L 
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IV  V 

T3a  ^  l  Ri  L  V  ^  R]_  “ 2 
L 


T3b 


2  V 

T3b  ^  R  T3a  -» 

1  1  L 


T4 


1  IV 

x4  „  fVL  „  A- 

L  L 


T5a 


„  1  V  V  IV 

T5a  %  Ro  L  L  ^  Ro  2 

L. 
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substantiated  that  the  above  terms  dominate  on  the  large  scale  (see, 
e.g.,  Haltiner,  1971). 

2.2.3  Determination  of  the  mesoscale  balance  equation 

As  in  Section  2.1.4,  we  use  for  mid-latitudes  and  the  mesoscale, 

L  'v  100  km,  Ro  m  1,  L/a  m  0.01,  H  m  10  km,  and  m  0.001.  With  these 
values,  the  relative  order  of  magnitude  of  the  terms  in  (2.25)  is 
indicated  in  Roman  numerals  underneath  each  term.  As  with  (2.19), 
the  difference  between  I  and  II  is  an  order  of  magnitude.  Because 
we  are  interested  in  a  divergent  initialization,  we  include  the  terms 
containing  the  divergent  component  and.  because  we  are  interested  in 
accuracy  to  second  order,  we  neglect  terms  labeled  III  (T2d,  T3b, 

T5c,  and  T6) .  Therefore,  the  final  mesoscale  balance  equation  is 
composed  of  terms  T2a,  T2b,  T2c,  T3a,  T4,  T5a,  and  T5b,  and  will  be 
referred  to  as  (2.26).  The  difference  between  the  synoptic  and 
mesoscale  balance  equations  is  the  relative  importance  of  the  Coriolis 
force  on  these  scales.  On  the  synoptic  scale,  the  beta  term  is  a 
first-order  term.  On  the  mesoscale,  the  beta  term  becomes  two  orders 
of  magnitude  smaller  than  the  first-order  terms. 

In  order  to  achieve  maximum  consistency  with  the  forecast  model, 
which  should  help  the  model  "remember"  the  divergent  component,  the 
balance  equation  will  be  applied  on  sigma  surfaces.  Sigma  is  defined 
by 


P  -  Pfc  P  -  Pt 


a  s 


P* 


(2.27) 
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where  p  is  the  pressure  of  the  sigma  level,  p  is  the  pressure  at  the 
model  top,  pg  is  the  surface  pressure,  and  p^  =  pg  -  Pt*  Now,  for 
Pt  =  0,  we  can  show 


3a  2  3Z 


With  (2.28),  term  T3a  of  (2.26)  becomes 
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and  T4  becomes 
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Therefore,  the  mesoscale  balance  equation  in  sigma  coordinates  is 
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3a  p. 
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where  V,  and  V  are  defined  on  pressure  surfaces.  That  is,  if 

<l>  X 

V  • V  »  0,  then  V  »V,  is  not  necessarily  zero  because  of  the  slope  of 

p  ~i|j  ’  a  -i(j 

the  sigma  surface.  In  this  thesis,  the  terms  nondivergent  and  divergent 
will  always  apply  to  pressure  and  not  sigma  surfaces. 
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2.2.4  Effect  of  neglecting  the  local  rate  of  change  of  divergence 

The  local  rate  of  change  of  divergence  was  neglected  in  the 
determination  of  the  balance  equation  in  Section  2.2.1.  This 
assumption  is  often  made  on  the  synoptic  scale.  The  scale  analysis  in 
Table  3  shows  that  T1  of  (2.25)  is  an  order  of  magnitude  smaller  than 
the  largest  terms.  We  will  show  that  this  assumption  is  reasonable 
on  the  upper  meso-B  scale. 

We  can  examine  the  effect  of  the  local  rate  of  change  of 
divergence  term  in  the  diagnosis  of  geopotential.  Let  the  entire 
divergence  equation  be 

V2*  =  -  .  (2.32) 

A  one-dimensional  analysis  is  sufficient  and  hence  we  can  write  (2.32) 
as 


3D 

3t 


BD 

Let  —  be  expressed  with  a  spectral  representation  as 
a  t 


(2.33) 


i<x 


(2.34) 


where  D  has  dimensions  of  s  and  k 


2tt/L  is  the  wavenumber.  Therefore, 


we  can  write  (2.30)  as 
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3Zij>  A  i<x 

— ir  =  -  D  e 


3x 


(2.35) 


Integration  of  (2.35)  twice  with  respect  to  x  yields 


D  i*x 

— j  e  +  c .  x  +  c„ 


(2.36) 


where  c^  and  C£  are  arbitrary  constants.  Therefore,  the  error  in 
geopotential  associated  with  the  neglect  of  this  term  is  on  the  order 
of 


<t>  ^ 

error  y2 


(2.37) 


Choose  a  typical  mesoscale  divergence  of  10  ^  s  The  mesoscale 

length  scale  L  ■v  10^  m  and  the  velocity  scale  V  'v  10  m  s  ^  imply  a 
4 

time  scale  of  10  s.  Therefore, 


D  'v  10 


-2 

s 


Now  (2.37)  becomes 


error 


With  (2.38)  and  L  'v  10  m,  (2.39)  gives 


^  0.25  m 


2 


-2 

s 


(2.38) 


(2.39) 


error 


(2.40) 
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Hence  for  the  local  rate  of  change  of  divergence  term  to  be  important 
on  the  upper  meso-B  scale,  large  values  of  divergence  are  required. 
Fankhauser  (1974)  reported  that  -rr  was  large  in  the  vicinity  of 

a  L 

a  squall  line  but  was  small  away  from  the  active  convection.  With 

3D 

this  result  and  (2.40),  we  conclude  that  neglecting  the  - —  term  on 

3 1 

the  mesoscale  is  an  acceptable  approximation.  The  term  might  become 
important  on  smaller  scales  and  may  have  to  be  included  there. 

2.3  Boundary  conditions  required  by  limited  domains 

The  diagnostic  equations  for  the  mesoscale,  derived  in  this 
chapter,  are  second-order  elliptic  PDEs.  For  purposes  of  this 
discussion  of  boundary  conditions,  they  can  be  represented  as  Poisson 
equations.  For  global  or  hemispheric  models,  boundary  conditions  have 
not  been  a  major  problem.  However,  the  attainable  resolution  has 
been  coarse  at  best.  Limited-area  models  have  been  developed  primarily 
to  obtain  increased  horizontal  as  well  as  vertical  resolution.  But, 
because  the  limited-area  domains  have  boundaries,  boundary  conditions 
must  be  specified  on  the  dependent  variables  or  else  the  problem  of 
solving  the  equation  is  not  mathematically  well  posed  (Ames,  1977). 

The  specification  of  boundary  conditions  for  limited-area  NWP  models 
and  their  initialization  schemes  has  indeed  been  a  problem.  In  this 
section  we  will  examine  the  boundary  conditions  required  by  the 
divergent  initialization  scheme  for  the  following  variables:  stream 
function,  geopotential,  omega,  and  velocity  potential. 
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2.3.1  Boundary  conditions  on  the  stream  function 

Several  of  the  terms  in  the  omega  and  divergence  equations 
require  knowledge  of  the  nondivergent  wind.  The  nondivergent  wind 
is  related  to  the  stream  function  by  (2.3a)  and  the  stream  function  is 
related  to  the  vorticity  field  by 

72*  =  t,  .  (2.41) 


Vorticity  is  calculated  at  the  interior  grid  points  and  ip  can  be 
calculated  on  the  interior  if  we  know  ip  .  There  are  numerous  methods 

D 

presented  in  the  literature  for  the  determination  of  '!>  (see,  e.g., 

15 

Phillips,  1958;  Anthes,  1976;  Brown  and  Neilon,  1961;  Bedient  and 
Vederman,  1964;  Hawkins  and  Rosenthal,  1965;  Sangster,  1960;  Shukla 
and  Saha,  1974;  Schaeffer  and  Doswell,  1979;  Stephens  and  Johnson, 
1978;  Endlich,  1967).  We  will  present  an  appropriate  method  after 
demonstrating  what  property  the  chosen  method  must  possess.  Before 
presenting  that  method,  however,  we  will  look  at  another  means  of 
avoiding  boundary  effects. 

A  method  which  may  improve  upon  the  results  obtained  by  the 
aforementioned  authors  has  been  used  by  Anthes  (1976),  Anthes  (1978), 
Keyser  (1978),  and  Elsberry  and  Ley  (1976),  among  others.  The  idea 
is  to  initialize  on  a  domain  larger  than  that  over  which  the  model 
will  produce  the  forecast,  thus  minimizing  boundary  effects  by  moving 
them  away  from  the  forecast  domain.  However,  we  choose  not  to  use 
this  approach  in  this  thesis  for  the  following  reasons: 
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(1)  It  is  not  always  possible  to  enlarge  the  domain. 

For  example,  in  certain  situations,  auxiliary  rawinsonde 
measurements  are  taken  over  limited  areas  (Hill  and 
Turner,  1977).  Hence,  there  are  no  data  on  which  the 
enlarged  initialization  domain  could  depend.  Also, 
it  is  not  always  desirable  to  enlarge  the  initialization 
domain.  Enlarging  the  domain  could  mean  incorporating 
ocean  or  other  data-sparse  regions  where  the  boundary 
values  of  the  meteorological  variables  contain  more 
uncertainty  than  in  data-dense  regions.  Thus  expanding 
the  initialization  domain  could  result  in  an  inferior 
initialization  on  the  forecast  domain. 

(2)  It  is  computationally  inconvenient  to  enlarge  'he  domain. 

(3^  We  will  show  that,  by  taking  a  few  reasonable 

precautions,  the  boundary  values  of  (and  x)  can  be 
determined  to  sufficient  accuracy  for  purposes  of  this 
thesis  without  expanding  the  initialization  domain. 

Keyser  (1978)  analyzed  the  effect  of  boundary  conditions  on 
the  solution  of  a  Poisson  equation  such  as 

V2<f(x,y)  =  F(x,y)  (2.42) 

where  the  rectangular  domain  is  defined  in  x  and  y  by  0  <  x  < 

and  0  y  <  D^.  Following  Morse  and  Feshbach  (1953),  that  portion  of 

the  solution  due  to  D  for  the  right  boundary  (q  )  can  be  written  as 

d  RB 
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.  t  mr 

smh  —  x  n 

«  m  D-  (D_ 

/  N  2  ®  1  mr  1  ,  /T  x  .  nn 

>RB(x’y)  =  d:  “7~~T“  sin  57  y  *(Li,y)  sin  57  y  dy 

1  n=l  smh  —  L1  1  Jo  1 

1  (2.43) 


There  are  four  other  terms  in  the  total  solution,  one  for  each  of  the 
other  three  boundaries  and  one  for  the  forcing  function.  For  this 
discussion,  we  need  to  consider  only  one  boundary.  Since  (2.42)  is 
linear,  we  may  interpret  (2.43)  for  a  single  harmonic  without  loss  of 


generality.  Define  y  as 


D1  L1 

sinh  L 

°1  1 


(2.44) 


For  a  given  distance  from  the  boundary,  for  increasing  wavenumber,  the 


value  of  y  decreases.  Therefore,  y  is  called  a  damping  factor.  Keyser 
plotted 


c,  1  x  ^ 

u  =  f(_ -  ,  — ) 

U1  L1 


(2.45) 


and  this  is  reproduced  as  Fig.  1.  From  this  argument  emerged  two 
very  important  conclusions  (Keyser,  1978): 


(1)  The  influence  of  boundary  conditions  on  the  solution 
decreases  exponentially  with  distance  from  the 
boundaries . 

(2)  It  is  most  important  to  accurately  specify  the  large- 


scale  (low-wavenumber)  variation  of  the  boundary 


Plot  of  contours  of  p  representing  the 
fractional  damping  of  the  boundary 
condition  at  x/L  =  1  as  a  function  of 
distance  normal  to  the  boundary  and  wave- 
number  (Keyser,  1978). 
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conditions  since  amplitude  errors  for  low-wavenumber 
components  of  the  boundary  conditions  damp  less  rapidly 
with  distance  away  from  the  boundary. 

Keyser  extended  the  analysis  to  a  discrete  domain  but  the  conclusion 
regarding  correct  specification  of  the  large-scale  variation  of  the 
boundary  conditions  was  unaltered. 

Therefore,  in  examining  the  adequacy  of  the  methods  for 
determining  ’J>g,  we  need  only  concern  ourselves  with  choosing  a  method 
capable  of  correctly  specifying  the  low- wave number  variation  of  \p  . 

D 

WAM,  Anthes  (1976),  and  Keyser  (1978)  report  that  when  Anthes  (1976) 

method  was  used  to  evaluate  <p  ,  the  required  correction  to  v 

B  n  , 

obs 

(observed  boundary  normal  wind  component)  was  on  the  order  of  a 

few  tenths  of  a  meter  per  second.  Since  the  magnitude  of  v^  itself 

was  about  two  orders  of  magnitude  larger  than  the  correction  applied, 

then  the  large-scale  features  are  retained  in  the  corrected  v  .  In 

n 

other  words,  since  the  observed  wind  was  predominantly  nondivergent , 

v  was  almost  v  .  We  conclude  that  Anthes  (1976)  method  may  be 
obs  ip 

used  to  solve  (2.41)  in  the  determination  of  the  nondivergent  wind 
for  the  scales  studied  here. 

Anthes  (1976)  method  is  based  on: 


dip 

9s 


(Vn 

obs 


+  CB) 


(2.46) 


where 
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c 


B 


r  t  vnK 

p  '  obs 


d 

5 


(2.47) 


and  where  L  is  the  length  of  the  domain  perimeter.  First,  c„  is 
P  B 

computed  from  (2.47).  Then,  after  one  value  is  specified,  (2.46) 

15 

is  integrated  around  the  boundary  to  obtain  a  complete  set  of  i>_  values. 

D 

That  is,  the  mean  divergence  over  the  domain  is  removed  from  v  by 

nobs 

applying  an  equal  correction  to  each  boundary  observed  normal  wind 
component . 


2.3.2  Boundary  conditions  on  geopotential  for  the  balance  equation 


We  know  from  the  discussion  of  the  previous  section  that  the 
large-scale  variation  in  geopotential  must  be  accurately  specified. 

The  usual  method  of  calculating  4>_  has  been  to  assume  that  the  boundary 

D 

winds  are  geostrophic  and  integrate 


3<t>  _  —  3tp 
3s  ~  t  3s 


(2.48) 


around  the  boundary.  WAM  reported  that  this  approach  was  unsatisfactory 
because  wind  analysis  errors  accumulated  around  the  boundary.  Therefore, 
WAM  used  the  observed  heights  on  the  boundary.  This  did  indeed 
preserve  the  large-scale  geopotential  variation. 

Keyser  (1978)  and  Anthes  (1978)  used  FNWC  (Fleet  Numerical 
Weather  Central,  Monterey,  California)  analyses  for  numerous 
initializations  and  subsequent  forecasts.  They  claimed  that  the  large- 
scale  geopotential  boundary  variation  was  preserved  and  hence 


4>_  =  4>„  was  an  adequate  boundary  condition.  That  is,  it  did  not 

B  B  , 

obs 

introduce  any  significant  error  in  geopotential  on  the  domain  interior 
Therefore,  we  adopt  the  use  of  observed  boundary  heights  as  the 
boundary  condition  for  the  balance  equation  in  this  thesis. 

2.3.3  Boundary  conditions  on  omega  for  the  omega  equation 

The  omega  equation  will  be  solved  by  three-dimensional  relaxation 
and  therefore  boundary  values  of  omega  for  the  top,  sides,  and  bottom 
of  the  domain  volume  are  required. 

For  the  top  of  the  domain  (250  mb),  we  set  omega  equal  to  zero. 

If  the  250-mb  level  is  the  top  of  the  features  and  circulations 
being  modeled,  then  this  choice  makes  meteorological  sense.  Since 
the  200-mb  level  is  the  top  data  level  used,  then  there  is  really  no 
other  reasonable  alternative.  Also,  the  model  itself  requires  zero 
omega  at  the  top  pressure  level. 

Omega  is  also  set  to  zero  on  the  side  faces  of  the  domain  volume. 
This  is  mathematically  expedient  but  probably  somewhat  unrealistic. 

The  meteorologically  interesting  feature  is  normally  placed  near  the 
center  of  the  domain.  Hence  there  are  usually  no  large  divergence 
values  near  the  boundaries.  With  this  precaution,  omega  values  of 
zero  at  the  side  boundaries  approximately  represent  the  correct  large- 
scale  boundary  variation  of  omega.  We  will  adopt  this  precaution  here 
The  PSU  model  uses  terrain  heights  as  the  lower  boundary 
condition  on  geopotential.  The  terrain  in  turn  induces  an  omega  at 
the  surface.  The  terrain-induced  omega,  wt,  can  be  written 
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=  _  SR 


RT  -i> 


V,-Vh 


(2.49) 


where  h  is  the  terrain  height.  We  use  scale  analysis  to  estimate  the 
importance  of  this  effect.  From  (2.11),  we  can  write 


Therefore,  (2.49)  becomes 


(2.50) 


=  -S—  V  -vh  —  —  =  — 

RT  -tj;  H  L  HL 


(2.51) 


Using  the  scales  H  %  10  km,  L  'v  100  km,  V  'v  10  m  s  \  and  a  terrain 

gradient  of  0.001  (a  100  m  rise  in  100  km  is  not  uncommon),  we  get 
•  -6  -1 

Zt  'v*  10  s  Note  that  in  mountainous  terrain,  a  slope  of  0.01 
would  be  more  appropriate  and  Z would  be  an  order  of  magnitude  larger. 

These  terrain-induced  vertical  velocities  are  indeed  significant. 
In  more  familiar  coordinates,  we  get 


•  2  -6  -1 

w  =  -  p  Z  -^10  cb  10  s 
t  st 


10-4  cb  s-1  -v  10  cb  d"1 


(2.52) 


Vertical  velocities  in  the  mid-troposphere  normally  don't  exceed 
10  cb  d  ^  except  in  sharp  troughs  and  precipitation  areas.  We  conclude 
that  the  terrain  effect  omega  values  should  be  used  as  the  lower 
boundary  condition  for  the  omega  equation  relaxation. 
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2.3.4  Determination  of  appropriate  boundary  conditions  on  velocity 
potential 

Velocity  potentials  are  required  so  that  the  divergent  wind 
component  can  be  supplied  to  the  omega  and  balance  equations.  Using 
(2.3b)  and  (2.13),  we  obtain 

V2X  =  D  (2.53) 


or 


»2X  -  -  If  •  (2.54) 

We  can  solve  (2.53)  given  values  of  omega  and  Xg- 

Several  methods  for  the  determination  of  Xg  have  been  used 
previously  (see,  e.g. ,  Brown  and  Neilon,  1961;  Bedient  and  Vederman, 

1964;  Shukla  and  Saha,  1974;  Endlich,  1967;  Schaeffer  and  Doswell, 

1979;  Stephens  and  Johnson,  1978).  Most  of  the  above  methods  involve 
at  least  one  assumption.  Two  of  the  methods  (Schaeffer  and  Doswell, 

1979;  Stephens  and  Johnson,  1978)  have  accurately  separated  observed 
wind  fields  into  nondivergent  and  divergent  components.  However, 
here  we  will  diagnose  vertical  velocities  and  therefore  divergence 
values.  We  want  a  divergent  component  corresponding  to  these  divergence 
values  and  not  necessarily  related  to  the  observed  or  nondivergent 
winds.  We  will  develop  a  means  of  accurately  specifying  the  large- 
scale  boundary  variation  of  velocity  potential. 


The  divergence  theorem  may  be  written 


V*G  dA  =  ®  n*G  ds 


(2.55) 


where  G  is  an  arbitrary  horizontal  vector,  A  is  the  domain  area,  and  s 
is  distance  along  the  domain  perimeter,  the  line  integral  being 
positive  in  the  counterclockwise  sense.  Applied  to  the  velocity 
vector,  (2.55)  becomes 


' 

' 

' 

’ 

(  - 

V-V  dA  = 

4 

D  dA  = 

>  n*V  ds 

• 

(2.56) 

In  words,  (2.56)  states  that  the  integral  of  the  divergence  over  the 
entire  domain  can  be  calculated  by  integrating  the  normal  boundary 
component  around  the  boundary.  That  is,  if  the  normal  velocity 
component  integrates  to  zero,  there  is  no  net  divergence  over  the 
domain. 

In  finite  difference  form,  (2.56)  becomes 


M-l  N-l  NBP 

Ax  I  l  D,  .  =  Z  v 

i=2  j=2  i,j  k-1  nk 


(2.57) 


where  M  is  the  number  of  grid  points  in  the  y  direction,  N  is  the 
number  of  points  in  the  x  direction,  and  NBP  is  the  total  number  of 
boundary  points.  Now,  from  the  omega  equation,  we  can  calculate  the 
left  hand  side  of  (2.57).  Therefore,  we  know  the  mean  value  of  the 
normal  velocity  component.  Since  ■  0,  it  follows  that 


64 


n*V  *  n-^X  =  f*. 
~  X  3n 


Now  we  know  the  mean  boundary  value  of  since 


(2.58) 


v  *  & 
n  3n 


(2.59) 


where  the  subscript  c  means  the  correct  value. 

The  complete  method  for  obtaining  x0  with  correct  low-wavenumber 

D 

variation  is  as  follows: 


(1)  Use  (2.57)  and  (2.59)  to  compute  the  exact  mean  value 
of  the  normal  derivative  of  velocity  potential  required 
to  satisfy  the  known  forcing  function. 

(2)  Set  x  over  the  domain  and  xB  equal  to  zero. 

(3)  Begin  to  solve 


7  x  =*  D 


by  relaxation  where 


72x  s  Xi,j+1  +  Xi,j-1  +  xi+l,J  +  Xi-l,j  I  A_Xi,j 

4 (Ax) 2 


where  Ax  is  the  grid  increment.  Apply  n  iterations. 

(4)  With  a  one-sided  difference,  compute  the  mean  boundary 

value  of  that  currently  exists  on  the  domain  after  n 
dn 


(2.60) 


(2.61) 


iterations 

3n 


,  where  the  subscript  a  means  the 


actual  mean  boundary  value). 


t 
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(5)  Compute 


BADJ  =  v 

n 


_  IX 

3n 


a 


where  the  LHS  of  (2.62)  is  the  mean  correction  that 

must  be  applied  to  the  existing  so  that  the  known 

dn 

boundary  value,  v  ,  will  be  realized. 

(6)  Apply  the  mean  normal  derivative  correction  by 

extrapolating  outward  from  the  first  interior  grid 
point.  That  is,  for  each  boundary  point,  compute 


XB  *  Xg.j  +  Ax- BADJ 

where  xB-1  is  the  value  of  x  at  the  first  interior 
grid  point. 

(7)  Return  to  step  (3).  Repeat  this  cycle  m  times. 

(8)  Xg  now  has  the  correct  large-scale  boundary  variation. 
Apply  a  direct  solver  (Rosmond  and  Faulkner,  1976; 
Swarztrauber  and  Sweet,  1975)  to  obtain  x  over  the 
entire  domain  for  the  given  x,>  values  and  forcing 

D 

function. 

The  advantages  of  this  procedure  are: 

(1)  It  imposes  no  arbitrary  boundary  conditions. 

(2)  The  method  is  physically  realistic  and  mathematically 
sound.  It  is  physically  realistic  because  information 


(2.62) 


(2.63) 
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about  the  known  forcing  function  is  transmitted 
to  xD  through  the  relaxation  procedure.  It  is 

D 

well-founded  mathematically  because  it  forces  the 
boundary  value  of  to  approach  the  value  it  must 
have  for  the  given  forcing  field.  We  are  interested 
in  the  correct  gradient  of  velocity  potential  and  not 
X  itself  since  only  the  x  gradient  has  physical  meaning. 

The  boundary  values  of  x  that  result  are  those  that 
would  exist  if  the  domain  were  infinite.  The  method 
is  mathematically  equivalent  to  using  a  Green's 
function  solution  of  (2.53)  (Hayek,  personal  communication; 
Morse  and  Feshbach,  1953).  The  applicable  Green's  function 
for  Cartesian  coordinates  is 

G(x,xo,y,yQ)  *  ^  In  [(x-xq)2  +  (y-yQ)2]  (2.64) 

The  solution  to  (2.53)  is  therefore  (Morse  and  Feshbach, 

1953) 

r  r 

X(x,y)  =  D(xQ,yo)  G(x,xo>y,yo)  dxQ  dyQ  (2.65) 

where  (x,y)  is  a  boundary  point  and  (x^.y^  is  any 
interior  point  where  the  divergence  D  is  defined.  For 
each  boundary  point,  (2.65)  must  be  applied  for  every 
interior  point.  Green's  function  solutions  have  the 
advantage  that  no  boundary  conditions  on  x  are  required. 
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However,  (2.64)  becomes  infinite  as  x  approaches  xq 
and  y  approaches  y  .  Therefore,  (2.65)  is  computationally 
difficult  to  use  in  practice.  Since  we  will  now 
demonstrate  that  the  method  presented  here  is  effective, 
we  will  use  it  in  this  thesis  for  the  determination  of 

XB* 

To  demonstrate  the  validity  of  the  Xt>  method,  several  numerical 

D 

experiments  were  performed  on  a  15  by  15  grid  with  a  grid  increment 

of  100  km.  The  scenario  was  to  generate  an  analytic  x  field,  compute 

the  divergence  (forcing  function)  at  each  grid  point,  and  calculate 

v  .  Then  the  x  values  over  the  entire  domain  were  zeroed  and  the 
n 

method  was  applied.  Table  4  contains  a  summary  of  the  experiments. 
Figures,  however,  are  included  only  for  Experiments  2.2a  through  2. 2d 
since  they  are  the  most  severe  meteorological  tests.  For  each  experiment 
the  value  of  n  in  step  (3)  was  nine.  That  is,  nine  relaxation  scans 
were  made  over  the  domain  for  each  boundary  normal  component  adjustment. 

In  Experiment  2.1,  a  circular  initial  velocity  potential  pattern 
was  placed  at  the  center  of  the  domain.  For  Experiment  2.2,  the 
center  of  the  circular  velocity  potential  pattern  was  shifted  left  to 
a  position  near  the  left  boundary  (Fig.  2).  For  Experiment  2.3,  the 
initial  x  field  was  given  by 
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Initial  velocity  potential  field  for 
Experiments  2.2a  through  2. 2d.  The 
contour  interval  is  1(P  s“^  corresponding 

to  a  divergent  wind  speed  of  1  m  s~  . 
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where  A  is  the  amplitude  of  the  sine  wave  and  here  corresponds  to  a 
divergent  wind  speed  of  5  m  s  \  Note  there  is  no  y  variation  of  x- 
In  Table  4,  MAXBDY  is  the  number  of  iterations  for  which  x„  was 

D 

corrected.  For  each  experiment.  Table  4  gives  the  v^  adjustment 
applied  (BADJ),  the  RMS  (root-mean-square)  error  in  u^  and  v  ,  the 
average  RMS  error,  the  percentage  reduction  in  RMS  error,  and  the 
number  of  iterations  (MAXBDY)  for  each  experiment  in  which  the  method 
was  applied.  The  percentage  reduction  in  RMS  error  is  the  percentage 
by  which  the  RMS  error  of  the  x  gradient  was  reduced  after  Experiments 
2.1a,  2.2a,  and  2.3a. 

In  Experiment  2.2a,  the  final  x  pattern  (Fig.  3a)  has  little 

resemblance  to  the  initial  pattern  (Fig.  2).  The  average  gradient 

error  corresponds  to  a  velocity  of  0.35  ns1.  In  Experiment  2.2b, 

with  the  method  applied  for  54  iterations,  the  final  x  pattern  (Fig.  3b) 

looks  much  closer  to  the  initial  pattern.  Experiments  2.2c  (Fig.  3c) 

and  2. 2d  (Fig.  3d)  show  additional  improvement  in  the  x„  accuracy  as 

o 

well  as  the  average  RMS  error. 

When  this  method  is  used  to  compute  Xg  and  subsequently  x>  the 
RMS  error  in  the  gradient  of  x  is  not  completely  eliminated.  Instead, 
the  error  decreases  and  then  oscillates  around  some  value  as  the  value 
of  the  boundary  correction  (BADJ)  levels  off.  Ior  example.  Experiment 
2.1b  actually  has  a  slightly  smaller  RMS  error  than  Experiment  2.1c. 

The  maximum  percentage  RMS  error  reduction  occurred  at  iteration  54  and 
oscillated  between  79  percent  and  84  percent  for  subsequent  iterations 
through  iteration  300.  Therefore,  it  seems  reasonable  to  establish  a 


I 
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c.  Experiment  2.2c.  d.  Experiment  2. 2d. 


Fig.  3.  Final  velocity  potential  field  for  Experiments  2.2a 
f  through  2. 2d.  The  contour  interval  is  10^  s--*- 

I  corresponding  to  a  divergent  wind  speed  of  1  m  s~l. 

r 

i 
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criterion  that  is  a  compromise  between  increased  accuracy  and  increased 
computational  time.  We  observe  that  when  BADJ  changes  by  less  than 
10  percent  over  its  previous  value  of  nine  iterations  earlier,  the 
majority  of  RMS  error  reduction  has  already  occurred.  BADJ  changed 
by  less  than  10  percent  at  iteration  number  45,  64,  and  64  for 
Experiments  2.1b,  2.2c,  and  2.3b,  respectively.  We  conclude  that,  when 
the  method  is  used  in  a  subsequent  chapter  to  determine  Xt>>  the 

D 

criterion  for  when  Xg  has  been  determined  to  sufficient  accuracy  is  when 

BADJ  changes  by  less  than  10  percent.  At  that  point,  a  direct  solver 

will  be  applied  to  obtain  the  velocity  potential  on  the  domain  interior. 

We  should  also  note  that  the  initial  x  patterns  for  Experiments  2.1 

(a,  b,  and  c)  and  2.2  (a  through  d)  are  severe  tests  of  the  method.  An 

average  v  of  almost  1  m  s  1  is  larger  than  values  normally 
obs 

encountered  (WAM;  Keyser,  1978). 

This  method  for  the  determination  of  Xt>  has  several  disadvantages: 

o 

(1)  The  method  is  an  iterative  technique  and  hence  requires 
more  computation  time  than  most  noniterative  techniques. 
However,  once  the  method  determines  a  set  of  xD  values 

D 

with  the  correct  low-wavenumber  variation,  iteration  is 
no  longer  required  and  a  direct  solver  is  used  to  obtain 
velocity  potential  on  the  domain  interior. 

(2)  Some  forcing  function  fields  require  more  normal  derivative 
boundary  corrections  (MAXBDY)  than  others.  This  is 
especially  true  if  small  values  of  the  forcing  function 
occur  near  the  boundary  or  if  a  large  range  of  forcing 


function  values  occur  over  the  domain  (both  conditions 
true  in  Experiments  2.2a  through  2. 2d).  We  normally 
do  not  place  a  boundary  in  the  vicinity  of  large  values 
of  the  forcing  function  (the  meteorologically  interesting 
feature).  That  precaution  will  help  minimize  the  effect 
of  this  disadvantage. 

To  summarize  this  section,  we  should  first  state  that  for  purposes 
of  divergent  initialization,  the  methods  which  have  been  used  previously 
may  be  unacceptable  because  they  do  not  accurately  specify  the  low- 
wavenumber  variation  in  x„.  We  have  developed  a  method  to  determine 

D 

Xg  in  such  a  way  as  to  insure  the  correct  large-scale  Xg  variation. 

As  previously  stated,  Keyser  (1978)  demonstrated  that  if  the  low- 
wavenumber  x0  variation  was  correct,  then  x  errors  due  to  x„  were 

D  D 

insignificant  in  the  domain  interior.  Therefore,  this  method  for  the 
determination  of  xD  in  conjunction  with  a  direct  solver  can  be  used  to 

D 

obtain  an  accurate  solution  of  the  Poisson  equation  for  velocity 
potential  on  a  limited  domain. 

2.4  Chapter  summary 

The  purpose  of  this  chapter  was  to  establish  the  divergent 
initialization  procedure  that  will  be  used.  Since  we  will  initialize 
the  PSU  model,  we  stated  which  forces  are  considered  in  the  model's 
equation  of  motion. 

In  the  first  section,  we  derived  a  complete  generalized  vertical 
velocity  equation.  Scale  analysis  was  applied  to  demonstrate  that  the 
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quasi-geostrophic  vertical  velocity  equation  resulted  for  the  synoptic 
scale-  We  then  established  the  vertical  velocity  equation  appropriate 
for  the  upper  meso-S  scale.  The  diabatic  term  was  an  important  term 
in  precipitation  areas  on  the  synoptic  scale  but  it  was  even  more 
important  on  the  upper  meso-S  scale.  In  fact,  on  the  upper  meso-B 
scale,  the  diabatic  term  dominated  the  other  forcing  functions.  Also, 
the  Coriolis  term  was  significantly  less  important  than  on  the  synoptic 
scale. 

In  the  second  section,  we  derived  the  complete  generalized 
divergence  equation.  The  nonlinear  balance  equation  was  shown  by 
scale  analysis  to  be  valid  on  the  large  scale.  The  balance  equation 
required  for  the  upper  meso-S  scale  was  established  and  presented  in 
sigma  coordinates.  On  the  upper  meso-S  scale,  the  Coriolis  term  was 
not  as  important  as  on  the  synoptic  scale.  We  showed  that  the  neglect 
of  typical  values  of  the  local  rate  of  change  of  divergence  was  an 
acceptable  assumption  on  the  upper  meso-B  and  larger  scales. 

In  the  last  section,  using  theoretical  and  mathematical  tools,  we 
established  what  boundary  conditions  should  be  used  for  stream  function, 
geopotential,  omega,  and  velocity  potential.  A  new  method  for  the 
correct  low-wavenumber  specification  of  Xg  was  presented  and  its 
effectiveness  was  demonstrated  for  severe  cases.  A  criterion  for  the 
termination  of  the  method  was  determined  experimentally.  This  new 
application  was  mathematically  equivalent  to  a  Green's  function 
solution  of  the  Poisson  equation. 
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3.0  THE  MESOSCALE  MODEL 


The  PSU  model  is  a  general,  predictive,  hydrostatic,  primitive 
equation,  meteorological  model  formulated  in  sigma  coordinates.  For 
a  complete  description  of  the  model,  see  AW  and  their  references. 

The  model  has  many  options  available  such  as  variable  terrain,  a 
moisture  cycle,  and  high-  and  low-resolution  PBL  physics.  The  model 
is  suitable  for  forecasting  flows  with  characteristic  horizontal 
wavelengths  of  about  10-2500  km  (meso-y  through  macro-B  scales) 
under  a  variety  of  meteorological  conditions.  It  is  indeed  a 
versatile  tool. 

3.1  General  description  of  the  model 

The  model  equations  described  in  AW  are  in  flux  form,  where  the 
vertical  coordinate  a  is  defined  by  (2.27).  A  Lambert  conformal  map 
projection  will  be  used.  There  are  equations  for  the  u  and  v  velocity 
components,  a  thermodynamic  equation,  and  continuity  equations  for 
mass  and  water  vapor. 

For  lateral  boundary  conditions  during  the  model  integration, 
a  linear  interpolation  in  time  between  the  balanced  conditions  at  the 
two  nearest  synoptic  times  is  used.  For  example,  for  a  00  GMT 
(Greenwich  Mean  Time)  to  06  GMT  forecast,  the  model  uses  boundary 
values  on  u,  v,  T,  $,  and  q  that  are  linearly  interpolated  in  time 
between  the  boundary  values  at  00  GMT  and  12  GMT.  These  "open" 


boundaries  allow  features  to  enter  the  domain  at  the  inflow  boundaries 


and  leave  the  domain  at  the  outflow  boundaries,  without  significant 
reflection  of  wave  energy. 

A  staggered  grid  is  employed  with  u  and  v  defined  at  points 

("dot  points")  midway  between  where  the  other  variables  are  defined 

("cross  points";  reference  Fig.  4).  At  the  lowest  sigma  level  (a  “  1) , 
•  • 

p^,  4>,  and  o  are  defined  while  p  and  a  are  specified  at  the  top 
(a  =  0)  sigma  level.  The  dependent  variables  themselves  (u,  v,  T, 
to,  <)>,  and  q)  are  defined  at  the  forecast  levels  whereas  at  the 
intermediate  sigma  levels,  a  and  the  vertical  fluxes  of  u,  v,  T,  and 
q  are  defined  (reference  Fig.  5). 

The  time-differencing  scheme  used  is  the  pressure-averaging 
technique  of  Brown  and  Campana  (1978) .  This  scheme  permits  a  larger 
time  step  while  meeting  the  linear  stability  (Courant-Friedrichs-Lewy 
or  CFL)  criterion  for  the  advection  term.  However,  in  preliminary 
forecasts  with  the  model,  time  splitting  (the  separation  of  the  odd 
time-step  solution  from  the  even  time-step  solution)  was  a  problem. 

As  a  result,  Anthes  (1978)  and  McNab  (unpublished)  incorporated  a 
low-pass  time  smoother  (Asselin,  1972;  Robert,  1966)  into  the  model. 
The  subsequent  model  performance  was  improved.  Unfortunately,  because 
of  numerical  stability  considerations,  the  time  smoother  requires  a 
smaller  time  step. 

The  model  employs  both  vertical  and  horizontal  diffusion.  In  the 
low-resolution  PBL  version  of  the  model  used  here,  a  simplified 
version  of  Deardorff's  (1972)  bulk  PBL  parameterization  is  used.  No 
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Fig.  4.  A  portion  of  the  staggered  horizontal  grid. 

The  horizontal  velocity  components  are 
defined  at  the  dot  points.  All  other  variables 
are  defined  at  the  cross  points. 


dr,  w'x‘ 


V,T,^,W,« 


Fig.  5.  Vertical  grid  structure  of  the  mesoscale  model 
showing  vertical  indexing  and  the  levels  at 
which  the  variables  are  defined. 
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other  levels  have  vertical  diffusion.  The  horizontal  diffusion  is 

t 

applied  to  each  variable  and  level  and  is  required  for  numerical 
stability  because  the  model  permits  nonlinear  interactions.  The 
horizontal  diffusion  scheme  is  that  used  by  Smagorinsky  et  al .  (1965) 


plus  a  constant  background  value  for  additional  smoothing  if  desired. 
Additionally,  Warner  and  McNab  (unpublished)  added  the  capability  to 
enhance  the  horizontal  diffusion  near  the  boundaries.  For  a  variable- 
size  elliptical  area  defining  the  domain  interior,  the  horizontal 
diffusion  is  not  enhanced.  Outside  that  area,  the  variable  portion 
of  the  diffusion  term  is  multiplied  by  a  coefficient.  The  value  of 
that  coefficient  is  one  at  the  edge  of  the  ellipse.  For  every  grid 
increment  of  distance  away  from  the  ellipse,  the  value  is  increased 
by  a  specified  amount  (variable  name  SPONGE).  That  is,  the  enhancement 
coefficient  is  equal  to  SPONGE  times  the  distance  in  grid  units  that 
the  grid  point  is  outside  the  ellipse.  Therefore,  the  horizontal 
diffusion  increases  as  the  boundary  is  approached.  However,  the  total 
value  of  the  horizontal  diffusion  coefficient  is  restricted  to  40 
percent  of  the  maximum  allowed  by  a  linear  stability  analysis. 

The  model's  moisture  and  cumulus  cloud  parameterizations  are  a 
simplified  version  of  Anthes'  (1977)  scheme. 

The  PSU  model  has  been  tested  and  verified  for  a  relatively  large 
number  of  cases  (WAM;  Anthes,  1978;  Shaginaw,  1979). 


3.2  Specific  model  parameters  used  in  this  thesis 

The  domain  chosen  is  a  30  by  35  grid  with  a  grid  increment  (Ax) 
of  120  km.  The  domain  is  centered  at  AIN  and  95W  and  therefore  covers 
most  of  the  contiguous  48  United  States  (US) .  The  time  step  required 
for  computational  stability  is  180  s. 

The  pressure  at  the  top  model  sigma  surface  is  fixed  at  250  mb. 
There  are  six  layers  between  the  seven  sigma  levels  of  0.0,  0.25,  0.4, 
0.55,  0.7,  0.85,  and  1.0. 

The  specific  value  of  the  surface  drag  coefficient,  C^,  to  be 
used  was  difficult  to  determine.  Anthes  (1978)  reported  that  the 
value  of  consistent  with  the  bulk  PBL  parameterization  should  be 
between  0.001  and  0.003.  Therefore,  a  value  of  0.002  will  be  used 
here  for  C^. 

The  variable  enhanced  horizontal  diffusion  scheme  was  tested  for 
SPONGE  values  of  0.0,  1.0,  5.0,  12.5,  and  25.0  for  a  120-km  Ax  and  the 
data  set  to  be  described  in  the  next  chapter.  In  the  version  used 
here,  the  size  of  the  ellipse  is  such  that  the  ellipse  passes  within 
4Ax  of  the  center  of  each  boundary.  The  value  of  SPONGE  chosen  was 
5.0.  That  is,  midway  between  the  corner  points  on  any  one  side  of 
the  grid,  the  enhancement  factor  would  be  20.  The  SPONGE  value  of 
5.0  was  chosen  because  it  represents  a  compromise  between  smoothing 
the  2Ax  noise  that  can  penetrate  the  domain  from  the  outflow  boundary 
and  not  noticeably  smoothing  the  large-scale  features  of  interest. 
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3.3  Static  initialization  procedure 


The  initialization  procedure  utilized  by  WAM  in  previous 
forecasts  will  be  briefly  covered  here.  Winds  were  subjectively 
analyzed  at  the  850- ,  700- ,  500- ,  400- ,  300- ,  and  200-mb  levels. 
Vorticity  fields  were  obtained  directly  from  these  wind  analyses 
using  (3.1). 


3 v  3 ii 

3 x  m  3y  m 


(3.1) 


In  (3.1),  x  and  y  are  the  horizontal  coordinates  of  the  Lambert 
conformal  map  projection  and  m  is  the  map  factor.  The  stream 
function  was  then  obtained  using 


(3.2) 


given  ,  the  stream  function  on  the  lateral  boundaries. 

B 

Geopotentials  were  calculated  from  the  nonlinear  balance  equation 


V2<t>  =  fv2ip  -  2m2  0l>xy  -  +  8l^y  +  (3'3) 

where  8,  =  and  y,  =  .  Of  course  to  solve  (3.3),  is  required. 

1  9y  I  3x  a 

The  observed  boundary  height  values  were  used  for  4>g  in  (3.3). 

Finally,  temperatures  were  calculated  from  the  derived  geopotentials 


through  the  hydrostatic  equation.  The  nondivergent  winds  and  balanced 
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temperatures  were  vertically  interpolated  from  the  levels  at  which 
they  were  calculated  to  the  model  sigma  levels. 
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4.0  SYNOPTIC  CASE  CHOSEN  FOR  STUDY 


We  will  use  data  for  two  synoptic  times.  For  each  time,  the 
available  data  consist  of  vector  winds  on  the  synoptic  rawinsonde 
network.  Shaginaw  (1979)  subjectively  analyzed  the  data  for  each  of 
the  six  standard  upper-level  pressure  surfaces  and  for  sea  level 
pressure.  He  then  manually  digitized  the  data  for  several  hundred 
points  at  each  level.  The  data  were  objectively  analyzed  with  a 
Cressman  (1959)  scan.  The  data  were  digitized  at  enough  points  so 
that  the  Cressman  scan  reproduced  well  the  sharp  gradients  and  smaller 
features. 

Shaginaw  (1979)  provided  a  comprehensive  discussion  of  the 
synoptic  situation  for  the  period  17-21  November,  1975,  and  that 
discussion  will  not  be  repeated  here.  However,  we  will  briefly 
summarize  the  meteorological  conditions  at  the  two  synoptic  times  used 
here.  Those  times  are  12191175  and  00201175,  where  the  hhddmmyy  format 
is  broken  down  in  hh  =  hour,  dd  *  day,  mm  =  month,  and  yy  =  year. 

For  example,  12191175  represents  12  GMT,  19  November  1975. 

A  high  pressure  area  (ridge)  persisted  over  the  southeastern  US 
between  the  times  12191175  and  00201175.  This  ridge  contributed  to 
a  severe  pollution  episode  at  Pittsburgh,  PA.  The  ridge  also  provided 
a  continuous  supply  of  low-level  moisture  from  the  Gulf  of  Mexico  to 
the  Great  Plains  and  adjacent  states.  This  moisture  supply  was  a  key 
factor  in  the  precipitation  which  occurred  between  12191175  and  00201175. 
The  moisture  initialization  used  was  provided  by  Wolcott  (1979) .  He 
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developed  a  scheme  incorporating  satellite,  surface,  and  rawinsonde 
data  into  the  relative  humidity  analysis.  The  initial  relative 
humidity  fields  for  sigma  levels  3%  and  5^  are  given  in  Fig.  6.  The 
satellite  picture  used  for  input  into  the  moisture  analysis  is  given 
in  Fig.  7.  It  is  important  to  note  that  the  initial  moisture  field 
is  saturated  over  a  large  portion  of  the  central  Plains  states  and 
the  upper  Midwest. 

4.1  The  synoptic  situation  at  12191175 

The  observed  sea  level  pressure  (SLP)  is  presented  in  Fig.  8a.  A 
ridge  dominates  the  eastern  United  States.  A  trough  extends  from  a 
1008-mb  low  near  Big  Bend  through  Minnesota. 

The  closed  low  at  500  mb  (Fig.  8b)  is  centered  at  the  "four 
corners"  region.  Although  the  ridge  over  the  Southeast  is  weakening, 
the  height  gradient  to  the  south  and  southeast  of  the  closed  low  is 
increasing.  Fig.  8c  is  the  observed  500-mb  wind  velocity.  There  are 
relative  minima  associated  with  the  closed  low  and  the  ridge  over 
the  Southeast.  Of  particular  interest  is  the  wind  maximum  entering 
the  domain  over  the  Calif ornia-Mexico  border.  This  wind  maximum 
("jet  streak")  will  proceed  to  move  around  the  south  end  of  the 
trough  and  up  the  east  side  resulting  in  significant  precipitation. 

Fig.  8d  is  the  observed  500-mb  temperature  field  at  12191175. 

Note  that  the  cold  tongue  extending  southward  from  Montana  reflects 


the  trough  position  well. 


Fig.  6.  Relative  humidity  analysis  valid  at  12191175  after  Wolcott  (1979).  Areas  of  greater  than  or 
equal  to  40,  80,  and  100  percent  relative  humidity  are  shown. 


Fig.  7.  Infrared-image  satellite  picture  of  the  domain  taken  at  1145  GMT,  19  November  1975. 
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Because  we  will  attempt  a  divergent  initialization  at  12191175, 
we  will  now  present  for  comparison,  the  results  of  the  nondivergent 
initialization  for  12191175  that  was  performed  as  outlined  in  Section 
3.3.  Fig.  9a  is  the  nondivergent  wind  speed  and  direction.  The 
differences  between  Figs.  9a  and  8c  in  the  vicinity  of  the  trough 
are  as  expected.  That  is,  the  wind  maximum  from  southern  California 
through  west  Texas  and  up  into  the  Plains  states  is  not  as  strong  in 
the  nondivergent  case.  Also,  the  horizontal  shear  over  New  England 
is  not  as  strong  in  the  nondivergent  case. 

Fig.  9b  is  the  balanced  height  field  obtained  from  the  non¬ 
divergent  winds  of  Fig.  9a.  When  compared  with  the  observed  heights 
(Fig.  8b),  the  balanced  heights  are  smoother,  the  trough  is  broader, 
the  height  at  the  center  of  the  closed  low  is  21  m  higher,  and  the 
height  gradient  southeast  of  the  trough  is  slightly  weaker. 

The  balanced  temperatures  are  given  in  Fig.  9c.  When  compared 
with  the  observed  temperatures  (Fig.  8d) ,  we  note  the  same  kinds  of  t 

differences  as  in  the  height  fields.  For  the  balanced  temperatures,  ; 

the  trough  is  warmer,  the  horizontal  temperature  gradient  north  of  I 

Mew  England  is  weaker,  and  the  ridge  over  the  Southeast  is  slightly  t 

i 

warmer. 

4.2  The  synoptic  situation  at  00201175 

Relatively  rapid  changes  occurred  between  12191175  and  00201175.  j 

A  low  developed  over  the  Texas  panhandle  and  moved  as  well  as  deepened  i 


rapidly  Co  a  1004-mb  low  over  eastern  Kansas  at  00201175  (Fig.  10a). 

The  trough  now  extends  from  south  Texas  through  the  Kansas  low  and 
then  through  Lake  Huron  and  into  New  England.  The  cold  air  has  now 
penetrated  into  the  northern  Rocky  Mountain  states  and  the  surface 
pressure  gradient  from  Wyoming  to  Kansas  is  much  stronger  than  12  hours 
before. 

From  the  observed  height  field  in  Fig.  10b,  we  see  the  closed 
low  has  moved  eastward  to  southwest  Kansas  at  a  speed  of  about 
17  ms  ■*".  The  height  of  the  low  center  fell  about  40  m.  The  height 
gradient  to  the  southeast  of  the  trough  has  strengthened.  This  is 
reflected  in  the  wind  speeds  and  directions  in  Fig.  10c.  The  wind 
maximum  previously  over  southern  California  has  increased  nearly 
12  m  s  1  while  moving  around  the  southern  end  of  the  trough  and  is 
now  centered  over  Abilene,  Texas.  Precipitation  occurred  on  the  cold 
side  and  ahead  of  this  jet  streak.  Also  note  that  the  wind  maximum 
north  of  Maine  has  increased  in  intensity.  The  position  of  this 
maximum  was  probably  responsible  for  the  precipitation  which  occurred 
in  the  western  Great  Lakes  region. 

The  temperature  field  (Fig.  lOd)  shows  an  increased  temperature 
gradient  to  the  south  and  southeast  of  the  trough  with  the  cold  air 
advancing  over  Kansas  and  Oklahoma  while  the  relatively  warm  air  still 
resides  over  the  Midwest. 


94 


4.3  The  observed  precipitation  amounts 

Precipitation  data  for  the  period  12191175  to  00201175  were 
extracted  from  the  National  Climatic  Center's  November,  1975, 

Hourly  Precipitation  Data  (Volume  25,  Number  11)  booklets  for  each 
state.  For  each  hour,  the  data  were  analyzed  at  grid  points  using  a 
Cressman  (1959)  scan  with  a  radius  of  influence  of  1.0. 

The  total  precipitation  for  the  first  three  hours  (12  GMT  to 
15  GMT,  19  November  1975)  is  given  in  Fig.  11a  and  t.  n  total 
precipitation  for  the  12-hour  period  is  given  in  Fig.  lib.  We 
observe  that  the  heaviest  12-hour  precipitation  (just  over  2.5  cm) 
occurred  from  south  central  Nebraska  to  northwestern  Kansas. 
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5.0  USE  OF  THE  OMEGA  EQUATION  IN  THE  DETERMINATION  OF  THE  DIVERGENT 
WIND  COMPONENT 


The  first  objective  of  this  chapter  is  to  review  previous  work 
on  the  use  of  the  omega  equation  to  diagnose  vertical  velocities. 

Then  we  will  derive  the  finite-difference  (FD)  version  of  (2.23),  the 
mesoscale  omega  equation.  The  FD  omega  equation  will  then  be  applied 
to  the  12191175  data  described  in  Chapter  4.  Experiments  will  be 
performed  to  confirm  the  validity  of  the  relative  order  of  magnitude 
of  the  terms  in  (2.23).  Finally,  sensitivity  experiments  will  be 
conducted  which  will  provide  information  on  how  accurately  the 
divergent  wind  component  can  be  determined. 


5.1  Previous  diagnostic  studies  using  the  omega  equation 


The  omega  equation  equivalent  to  that  of  Krishnamurti  (1968a) 
can  be  written: 


2 

V2 (c  w)  +  f2  ~  =  f  —  V  *75  +  f  ~  V  *V? 
s  3p2  3p  ~,p  Sp 


I, I, I  II, II, III 


+  -  V2V  •  VT  +  -  V2V  •  VT  -  f  |-  CD  +  f  u|5. 
P  P  -X  3p  9p  5p 


1,1, 1  II, II, III  II, II, II  II, II, III 


* 


A  %• 
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+  f  |r  70) -v  72q  +  f  i-  g  i_  v  x  T 

dp  dp  CpP  3p  9p 

T9  T10  Til 

II, II, III  I, I, I  II, II, II 


JL  v2h  _  2  A-  3_  J(l±  ii) 
s  ^  3t  3p  Jc3x  ’  gy 


CPP 
T12 

X,IV,X 


T13 

X,II,X 


o  i_  i_ 

3p  3y  3t 

T14 


X,III,X 


(5.1) 
(cont. ) 


where  H  is  the  sensible  heat  added  per  unit  time  from  a  water  surface, 
s 

Note  that  terms  T1  through  Til  correspond  exactly  to  those  of  (2.23) 
except  that  Krishnamurti  assumes  £<<f  in  T2.  Under  each  term  in 
(5.1),  the  first  Roman  numeral  is  the  order  of  magnitude  of  that 
term  from  the  scale  analysis  of  Chapter  2  (X  indicates  not  given). 
Krishnamurti  found  that  T12  was  small  except  over  large  water  bodies 
and  hence  T12  will  not  be  discussed  here.  Terms  T13  and  T14  cor¬ 
respond  to  terms  T13  and  T14,  respectively,  of  (2.19),  and  were  assumed 
to  cancel  in  (2.23)  by  virtue  of  the  geostrophic  vorticity  assumption. 

Krishnamurti  (1968b)  applied  (5.1)  on  a  2.5-degree  latitude 
by  2.0-degree  longitude  grid  for  four  synoptic  times  covering  the 
development  of  a  cyclone  in  mid-latitudes.  To  examine  the  magnitude 
of  omega,  he  chose  eight  points  in  the  vicinity  of  a  500-mb  trough 
at  one  of  the  four  synoptic  times.  First,  the  omega  values  at  the 
eight  points  were  computed  and  then  averaged.  From  this  average 
value,  the  relative  order  of  importance  to  the  total  omega  of  each 
term  is  given  as  the  second  Roman  numeral  under  each  forcing  function 
in  (5.1).  We  note  that  the  relative  importance  of  each  term  found 


by  Krishnamurti  is  exactly  that  determined  by  the  scale  analysis 
for  synoptic  scales  (Krishnamurti ' s  average  grid  increment  was  about 
200  km).  He  also  pointed  out  that  T8  and  T9  tend  to  cancel  and,  to 
a  lesser  extent,  so  do  T4  and  T6.  The  largest  vertical  velocity 
found  at  any  point  was  at  a  point  where  the  latent  heating  con¬ 
tribution  was  the  largest. 

Baumhefner  (1968)  used  Krishnamurti' s  (1968a)  diagnostic  model 
in  the  tropics  on  a  2.0-degree  latitude  by  2.0-degree  longitude  grid. 
He  studied  an  easterly  wave  spanning  four  synoptic  times  in  August, 
1961.  The  third  set  of  Roman  numerals  under  the  terms  on  the  RHS 
of  (5.1)  were  Baumhefner 's  results  at  the  500-mb  level.  We  see  that 
terms  T4,  T6,  T8,  T9,  and  T13  were  found  to  be  an  order  of  magnitude 
smaller  than  Krishnamurti’ s  results. 

Hawkins  (1972)  developed  a  diagnostic  model  for  computing  omega 
which  he  applied  to  three  disparate  cases  on  a  206-km  grid  mesh. 

He  reported  that  the  omega  values  obtained  from  the  complete, 
relatively  sophisticated  model  were  similar  to  the  first-guess 
values  diagnosed  with  the  quasi-geostrophic  omega  equation  with  the 
latent-heating  term  added.  That  is,  Hawkins'  findings  agreed  with 
the  other  researchers. 

In  summary,  we  conclude  that  previous  research  supported  the 
scale  analysis  used  in  Section  2.1  in  the  determination  of  the  omega 
equation  appropriate  on  the  synoptic  scale.  When  the  FD  form  of  the 
mesoscale  omega  equation  is  applied  to  the  12191175  data,  we  expect 
the  relative  importance  of  the  terms  will  be  the  same  as  that 
indicated  by  the  scale  analysis  for  the  mesoscale. 
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5.2  Derivation  of  the  finite-difference  (FD)  form  of  the  omega 
equation 


Before  we  can  present  or  discuss  FD  equations,  we  need  to  define 
some  FD  operations.  We  use  Shuman  and  Hovermale's  (1968)  notation 
for  the  FD  operators 


ax  S  (5'2a> 

*y  5  -  al-4,J>/ay  (5-2b) 

71  !  2  <5-2c) 

7  2  2  ai-l,,J>/2  (5-2d> 


where  j  is  the  east-west  index  and  i  is  the  north-south  index. 
We  will  use  the  "four-point"  operators 


“y  s  <Vi,j  +  2ai,j  + 


s  <al,Jrt  2  2“i,j  2  “i ,!-l>/4 


Vertical  differences  and  averages  are  defined  by 


H  (ak+4  +  V«J)/2 


6a  S  (ak4^  - 


(5.3a) 

(5.3b) 


(5.4a) 

(5.4b) 
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Omega  will  be  calculated  on  a  30  by  35  by  7  grid.  Each  30  by 
35  horizontal  section  has  a  grid  increment  of  120  km  and  is  centered 
at  41N  and  95W.  Fig.  12  illustrates  the  vertical  structure  of  the 
omega  equation  domain.  Note  that  velocity,  temperature,  omega,  and 
static  stability  are  defined  at  standard  constant-pressure  levels 
while  the  forcing  functions  for  the  three-dimensional  relaxation  are 
defined  at  the  mid-levels. 

Because  the  derivation  of  the  FD  form  of  the  omega  equation  is 
long  and  detailed,  it  is  presented  in  Appendix  1. 

5.3  Application  of  the  FD  omega  equation  to  the  12191175  data 

In  this  section  we  will  use  the  FD  omega  equation  to  compute 
omega  values  for  the  12191175  data  set.  We  will  first  examine  the 
diabatic  term.  Then  we  will  try  to  determine  to  what  accuracy  the 
vertical  velocities  can  be  determined. 

5.3.1  The  diabatic  term  and  the  parabolic  omega  profile 

As  described  in  Appendix  1,  we  use  a  parabolic  omega  profile 
and  the  observed  rain  rate  to  calculate  omega  values  due  to  the 
diabatic  term  only.  Hereafter,  these  omega  values  will  be  termed 
convective  omegas. 

Since  precipitation  amounts  were  measured  hourly  at  synoptic 
stations,  and  since  a  rainfall  rate  valid  at  12191175  was  desired, 
rainfall  amounts  from  2  hours  before  and  2  hours  after  12191175  were 
averaged  to  obtain  a  more  representative  rainfall  rate.  Fig.  13  presents 
that  precipitation  rate  in  cm  d 
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Fig.  12.  Vertical  structure  of  grid  used  to  calculate  the  omega 
field.  The  forcing  functions  F  are  calculated  at  half 
levels  and  omega  is  obtained  at  the  standard  pressure 
levels. 


Fig.  13.  Observed  rainfall  rate  in  cm  d  valid  at  12191175. 

This  rate  was  derived  from  the  observations  for  the 
period  10191175  to  14191175.  The  contour  interval 
is  1  cm  d-^. 
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We  should  mention  here  that  observed  precipitation  rates,  the 
input  data  for  the  omega  equation  diabatic  term,  can  be  obtained  from 
other  than  rainguage  measurements.  For  the  eastern  US,  rain  rates 
can  be  obtained  from  NWS  manually  digitized  radar  (MDR)  data  (Moore 
et  al. ,  1974).  Over  the  oceans,  meteorological  satellites  have 
provided  precipitation  rate  observations  (Adler  and  Rodgers,  1977).  A 
scanning  microwave  radiometer  on  the  satellite  measures  a  "brightness" 
temperature.  The  dominating  factor  in  the  determination  of  this 
temperature  over  water  is  liquid  water  drops  of  rainfall  size.  The 
brightness  temperatures  are  then  translated  into  rainfall  rates  using 
previously  derived  rainfall  rate-brightness  temperature  relationships. 
In  the  future,  satellites  may  also  provide  rainfall  rates  over  land. 
Therefore,  satellites  provide  a  means  of  obtaining  data  required  for 
the  diabatic  term. 

We  will  use  the  rainfall  rate  in  Fig.  13  in  the  solution  of  the 
finite-difference  form  of  the  omega  equation.  Before  solving  the 
entire  equation,  we  will  conduct  four  experiments  to  determine 
convective  omegas.  The  purpose  of  these  experiments  is  to  determine 
the  effect  of  uncertainty  in  the  observed  rainfall  rate  and  the  effect 
of  including  a  terrain- induced  omega  as  the  lower  boundary  condition 
for  the  parabolic  omega  profile. 

Table  5  gives  a  summary  of  the  convective  omega  experiments.  The 
column  labeled  "RAMT"  is  the  fraction  of  the  rainfall  amount  of  Fig.  13 
that  was  used  for  that  particular  experiment.  Fig.  14  is  the  500-mb 
convective  omega  field  for  Experiment  5.1.  The  purpose  of  Experiments 
5.1  through  5.3  was  to  determine  if  there  exists  a  significant 
variability  in  the  omega  values  calculated  for  the  computational  domain 


at  500  mb  when  the  precipitation  rate  was  altered.  Since  precipitation 
amounts  over  large  areas  are  determined  by  very  few  observations  in 
time  and  space,  the  precipitation  amounts  themselves  at  grid  points 
possess  a  sizable  uncertainty.  This  uncertainty  is  largest  in 
areas  of  strong  convective  activity.  A  25  percent  difference  between 
the  "observed"  precipitation  amount  and  the  precipitation  amount 
that  would  be  representative  of  a  grid  square  can  easily  be  imagined. 
For  a  25  percent  error  in  precipitation  rate  (represented  by 
Experiment  5.3),  there  was  approximately  a  20  percent  error  in  the 
500-mb  RMS  omega  value  over  the  entire  domain.  Therefore,  the  signi¬ 
ficance  of  the  diabatic  term  was  critically  dependent  on  the  accuracy 
of  the  rainfall  amounts  themselves.  Experiment  5.4  demonstrated  that 
when  terrain- induced  omega  values  were  not  used  as  the  lower  boundary 
condition,  there  was  only  a  small  effect  at  500  mb  but  there  was  a 
significant  effect  at  850  mb.  We  conclude  that  the  terrain  effect 
did  indeed  produce  significantly  different  omegas  in  the  lower  levels. 
Terrain- induced  omegas  should  be  used  since  they  are  consistent  with 
the  forecast  model  which  contains  variable  terrain. 

5.3.2  Experiments  with  various  terms  in  the  FD  omega  equation 

Experiment  5.5  consisted  of  obtaining  omega  values  using  the 
finite-difference  version  of  (2.23),  the  mesoscale  omega  equation. 

The  resultant  500-mb  omega  values  are  presented  in  Fig.  15.  From 
the  figure  we  can  see  the  strong  influence  of  the  diabatic  term  over 
eastern  Colorado,  the  Texas  panhandle  to  southwestern  Nebraska,  and 
over  Lake  Superior  and  vicinity. 
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Several  experiments  were  performed  to  support  the  scale  analysis 
used  In  obtaining  (2.23).  For  example,  experiments  were  completed 
with  the  quasi-geos  trophic  omega  equation,  with  differential  vorticity 
advection  only,  and  with  the  Laplacian  of  temperature  advection  only. 
From  these  experiments  the  following  conclusions  were  reached: 

(1)  The  omega  equation  experiment  without  a  diabatic  term 
greatly  underestimated  the  omega  values  in  the  precipitation 
areas.  In  fact,  the  diabatic  term  was  the  largest  single 
term  in  the  areas  of  precipitation.  Therefore,  a 
diabatic  term  must  be  included  on  the  mesoscale. 

(2)  The  quasi-geostrophic  omega  equation  with  a  diabatic 
term  overestimated  the  omega  values  in  the  precipitation 
areas. 

(3)  An  experiment  with  random  uncertainty  in  the  wind  field 
produced  only  small  RMS  changes  in  the  omega  values. 

(4)  The  scale  analysis  in  Section  2.1.4  was  supported.  That 
is,  the  various  forcing  functions  had  the  expected  relative 
influence  on  the  omega  values.  Therefore,  (2.23)  is  an 
appropriate  form  of  the  omega  equation  for  the  mesoscale. 

(5)  An  experiment  was  conducted  in  which  the  precipitation 
rate  used  for  the  diabatic  term  was  25  percent  less  than 
that  for  Experiment  5.5.  At  500  mb,  the  RMS  difference 
between  this  experiment  and  Experiment  5.5  was  15  percent 
of  the  mean  absolute  value  of  omega  for  the  entire  field. 

At  600  mb,  this  produced  a  10  and  15  percent  RMS 
difference  in  u^  and  v  ,  respectively.  That  is,  with  a 
reasonable  estimate  of  the  uncertainty  inherent  in 
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precipitation  observations,  the  percentage  change  in  the 
divergent  wind  components  was  almost  as  large. 

The  omega  values  used  to  determine  the  divergent  wind  component 
were  those  from  Experiment  5.5.  Then,  to  obtain  velocity  potential, 
the  method  described  in  Section  2.3.4  was  applied  with  the  boundary 
normal  derivative  adjusted  for  100  iterations.  The  resultant 
divergent  wind  fields  for  the  925-,  775-,  350-,  and  250-mb  levels 
are  given  in  Fig.  16.  In  Fig.  16a,  there  is  a  narrow  zone  of  con¬ 
vergence  from  central  Wisconsin  to  southwestern  Kansas  and  on  toward 
the  south.  This  same  pattern  is  evident  at  the  775-tnb  level  (Fig. 
16b) .  The  low-level  convergence  is  supported  by  the  divergence  aloft 
depicted  in  Figs.  16c  and  16d.  The  zone  of  divergence  extends  from 
Lake  Superior  to  western  Kansas  and  to  the  south.  The  divergent 
wind  fields  appear  to  be  vertically  consistent.  In  Chapter  7,  we 
will  see  if  the  region  of  maximum  vertical  motion  is  reflected  in  the 
forecast  started  from  a  divergent  initialization. 

In  the  next  chapter,  the  FD  form  of  the  balance  equation  which 
will  be  used  in  the  divergent  initialization  will  be  derived. 
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6.0  DERIVATION  OF  THE  FINITE-DIFFERENCE  (FD)  BALANCE  EQUATION 

The  purpose  of  this  chapter  is  to  derive  the  FD  form  of  the 
balance  equation  to  be  used  for  the  divergent  initialization  of  the 
PSU  model.  The  balance  equation  in  sigma  coordinates  appropriate 
for  this  purpose  is  given  by  (2.31). 

At  this  point  we  will  examine  how  balance  equations  in  general 
have  been  used  in  the  initialization  of  numerical  models.  The  common 
practice,  especially  on  the  large  scale,  has  been  to  use  observed 
geopotentials  and  solve  the  balance  equation  for  stream  function. 

For  example,  this  procedure  was  used  by  Sundqvist  (1975).  There 
are,  however,  at  least  two  reasons  why  an  alternative  approach  should 
be  employed. 

First,  Paegle  and  Paegle  (1976)  reported  that,  for  the  year 
beginning  1  June  1969,  virtually  every  Northern  Hemisphere  200-mb 
chart  exhibited  nonelliptic  geopotential  data  at  some  points.  The 
nonelliptic  geopotentials  were  especially  prevalent  in  the  summer 
above  severe  weather  events.  The  balance  equation  was  not  solvable 
in  regions  where  the  geopotentials  are  nonelliptic  and  some  approxi¬ 
mation  must  be  made  there. 

Second,  WAM,  Fankhauser  (1974),  and  AW  have  pointed  out  that 
typical  errors  in  surface  pressure  and  rawinsonde  temperatures  led 
to  large  errors  in  the  horizontal  gradients  of  geopotential  heights 
on  the  mesoscale  (10-100  km).  Fankhauser  reported  that  errors  in 
height  observations  of  more  than  30  m  were  common  on  the  mesoscale. 
Therefore,  we  will  follow  WAM  and  supply  the  nondivergent  wind  to  the 
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balance  equation  and  solve  for  geopotential.  We  now  proceed  to  the 
derivation  of  the  FD  balance  equation. 

One  possible  approach  for  the  FD  form  of  the  balance  equation 
would  be  to  write  each  term  of  (2.31)  directly  in  its  FD  form. 

Although  that  approach  is  conceptually  simple  and  relatively  easy 
to  do,  it  will  not  be  used  here.  One  objective  of  any  balancing 
scheme  should  be  to  retain  the  maximum  degree  of  consistency  with 
the  model  being  initialized.  This  should  help  minimize  noise  generated 
during  the  adjustment  phase  of  the  prediction.  That  is,  the  more 
consistency  between  the  initialization  scheme  and  the  model,  the 
nearer  the  balanced  initial  conditions  are  to  those  conditions 
exactly  compatible  with  the  model.  Therefore,  we  will  derive  the  FD 
form  of  the  balance  equation  directly  from  the  model  FD  equations 
themselves.  In  this  chapter,  we  will  present  the  PSU  model  FD 
equations,  derive  the  consistent  FD  balance  equation,  and  outline  how 
the  equation  should  be  used  in  the  initialization  procedure. 

6.1  The  PSU  model  FD  equations 

AW  provide  the  complete  set  of  the  FD  model  equations,  and 
most  will  not  be  repeated  here.  We  present  only  those  which  are 
essential  to  this  chapter. 

The  finite-difference  equations  associated  with  the  u  and  v 
component  equations  of  motion  (minus  those  terms  which  do  not 
correspond  to  a  term  in  (2.28))  are: 
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where  the  terms  p*u  and  p*v  represent  up*  and  vp*  '  ,  respectively. 
The  FD  form  of  the  continuity  equation  is 
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The  FD  equation  relating  omega  and  a  is 
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Finally,  the  hydrostatic  equation  is 


3  $ 

31n(a  +  p  /p.) 
t 


— CT 

RT 


(6.5) 


6.2  Derivation  of  the  FD  balance  equation  for  the  PSU  model 

We  know  that  a  divergence  equation  can  be  formed  from  the  u  and 

£ 

v  equations  of  motion  by  operating  on  the  u  equation  with  —  and  on 

the  v  equation  with  —  and  summing  the  result.  To  form  the  FD 

3y 

balance  equation,  the  procedure  must  be  slightly  modified.  In  the 
p^u  and  p^v  component  equations,  each  individual  term  is  defined 
at  a  dot  point  because  p^u  and  p*v  are  defined  at  dot  points.  However, 
we  will  use  the  balance  equation  to  solve  for  geopotential  and  on  the 
staggered  grid,  geopotential  is  defined  at  cross  points.  Therefore, 
in  deriving  the  balance  equation,  we  must  make  a  modification  such 
that  the  terms  are  defined  at  cross  points  rather  than  dot  points. 

The  result  of  differentiation  of  T1  in  (6.1)  with  respect  to  x 
is  Tl^  in  FD  notation.  However,  the  resultant  derivative  applies 
midway  between  dot  points.  In  order  for  the  derivative  to  apply  at 
the  desired  cross  point,  T1  must  first  be  averaged  in  the  y  direction 
(FD  operation  Tl^) .  For  T1  at  point  (i,j),  the  result  of  (Tl7)x  applies 
at  cross  point  (i+4,j+^)  as  desired,  where  integer  values  of  i  and  j 
refer  to  dot  points.  To  summarize,  the  desired  balance  equation  is 
derived  by  first  dividing  (6.1)  and  (6.2)  by  mp^.  Eqs.  (5. 2d)  and 
then  (5.2a)  are  applied  to  (6.1)  and  that  result  is  added  to  the 
result  of  applying  (5.2c)  and  then  (5.2b)  to  (6.2).  Transposing  the 


terms  containing  geopotential  to  the  LHS  of  the  equation  yields 
equation  (6,6),  the  complete  FD  balance  equation. 
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For  a  nondivergent  initialization  on  sigma  surfaces,  the  u  and  v 
appearing  in  terms  T3  through  T6  represent  the  nondivergent  wind 
components  u^  and  v^,  respectively.  For  a  divergent  initialization 
on  sigma  surfaces,  terms  T3  through  T6  are  each  applied  three 


because  the  u  and  v  in  those  terms  represent  the  total  wind  com,  '^er, t.s. 
That  is. 


u  *  u ,  +  u 

X 


(6.7a) 


V  =  V,  +  V 

<P  X 


(6.7b) 


For  example,  if  we  use  functional  notation  to  express  T3  as  a  function 
of  u  and  p^u,  then  we  can  write  T3  as 


T3  »  f(u,p.u) 


Substitution  of  (6.7a)  into  (6.8)  produces 


T3  -  ftu^.P*^)  +  f(u^,p*UY)  +  f(uv,p*Ujj) 


+  f(ux,Pll[ux) 


rr* x' 


x,r*  r 


The  fourth  term  in  (6.9)  is  neglected  because  the  scale  analysis  of 
Chapter  2  indicates  that  it  is  two  orders  of  magnitude  smaller  than 
the  first  term.  The  three  applications  of  T3  required  by  divergent 
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initialization  on  sigma  surfaces  are  represented  by  the  first  three 
terms  of  (6.9).  We  write  these  terms  as 

T31  *  f<vp*V 

T32  =  f(u^,p*ux) 

T33  -  f(ux,p*u^)  . 

Similarly,  for  terms  T4  through  T6,  we  can  write 
T41  -  f(«yp*vt> 

™  *  f<VP*V 


T62  «  f  (v,^,P*vx) 

T63  -  f(vx>P*vt)  • 

It  is  important  to  realize  that  the  "nondivergent"  or  "divergent" 


velocity  components  refer  to  pressure  surfaces  even  though  the  model 
is  formulated  on  sigma  surfaces.  We  calculate  the  velocity  components 
on  pressure  surfaces  and  interpolate  them  to  sigma  surfaces.  We  do 
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not  compute  x  °r  ♦  on  sigma  surfaces.  Thus  in  the  "nondivergent" 

initialization,  V  *v  *  0  but  7  -V  does  not  equal  zero  because  of  the 

P  -  o  - 

slope  of  sigma  surfaces  over  variable  terrain. 

Performing  the  above  FD  operations  on  the  terms  of  (6.6)  is  a 
long,  complex,  and  tedious  enterprise  and  is  therefore  presented  in 
Appendix  2. 


6.3  Application  of  the  FD  balance  equation 


To  apply  (6.6),  we  need  o  at  a  levels  and  cross  points,  T  at 

half-sigma  levels  and  cross  points,  and  u  ,  u  ,  v  ,  ,  v  ,  m,  and  f  at 

<P  X V  X 

half-sigma  levels  and  dot  points. 

The  available  data  consist  of  surface  observations  and  observations 
on  constant-pressure  surfaces.  Therefore,  the  data  will  be  inter¬ 
polated  to  sigma  surfaces  before  applying  (6.6).  This  will  provide 

T  at  half-sigma  levels.  We  obtain  u  and  v  at  half-sigma  levels  by 

X  X 

interpolating  the  divergent  components  computed  in  Chapter  5  from 
constant  pressure  to  sigma  levels.  We  will  use  the  observed  u  and  v 
at  pressure  levels  in 
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to  obtain  <p  at  the  mandatory  pressure  levels.  The 
divergent  wind  components  are  then  interpolated  to 
sigma  surfaces. 

3P* 

To  determine  o,  we  first  calculate  -  from  a 
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form  of  the  vertical  integral  of  mass  divergence.  Then,  for  each 
layer,  we  apply 


118 


•  ,  **  ,  — xy  x  .  — xy 

°k  -  (3T  +  “S»s  p*  +  V* 

x 


(6.11) 


to  obtain  a  at  the  full  intermediate  sigma  levels.  By  definition, 
a  is  2ero  at  the  top  and  bottom  sigma  levels.  We  have  now  defined 
all  of  the  quantities  necessary  to  use  (6.6). 

Once  geopotentials  have  been  obtained  at  all  the  half— sigma 
levels,  (6.5)  is  used  to  obtain  the  balanced  temperatures.  The 
geopotential  at  the  ground  is  known  from  the  terrain  elevation. 

The  balanced  temperatures  are  then  interpolated  to  the  half-sigma 
levels  and  the  initialization  process  is  completed. 


7.0  DIVERGENT  INITIALIZATION  APPLIED  TO  A  REAL  DATA  CASE 


The  purpose  of  this  chapter  is  to  discuss  in  detail,  five  12-hour 
forecasts,  all  of  which  began  at  12191175.  Forecast  1  is  initialized 
with  unbalanced  data.  The  observed  winds  and  temperatures  are 
interpolated  from  pressure  to  sigma  surfaces  and  are  inserted  into 
the  model  with  no  other  processing.  Forecast  2  is  initialized  with 
unbalanced  data  but  differs  from  Forecast  1  in  that  the  winds  are 
nondivergent .  The  nondivergent  winds  and  observed  temperatures  are 
interpolated  from  pressure  to  sigma  surfaces.  Forecast  3  is  initialized 
by  the  nondivergent  balance  method  described  in  Chapter  3.  That  is, 
geopotentials  are  calculated  on  pressure  surfaces  via  the  balance 
equation,  temperatures  are  derived' hydrostatically ,  and  these 
temperatures  are  interpolated  to  sigma  surfaces.  Forecast  4  follows 
a  nondivergent  initialization  on  sigma  surfaces.  Geopotentials  are 
calculated  on  sigma  surfaces  using  only  the  nondivergent  wind. 

Balanced  temperatures  are  derived  and  interpolated  to  sigma  surfaces. 
Forecast  5  follows  a  divergent  initialization  on  sigma  surfaces.  It 
is  similar  to  Forecast  4  except  that  the  terms  in  the  balance 
equation  containing  the  divergent  wind  components  are  also  used  and 
the  divergent  wind  component  obtained  in  Chapter  5  is  added  to  the 
nondivergent  wind  to  obtain  the  total  wind  field.  Forecasts  1  and  2 
are  considered  control  experiments  because  they  are  unbalanced. 

Forecasts  3,  4,  and  5  are  balanced  in  different  ways.  Figs.  17 
through  20  are  flow  charts  illustrating  how  the  various  initializations 
proceed.  Fig.  17  shows  how  the  nondivergent  wind  components  are 
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Fig.  17.  Flow  chart  depicting  the  calculation  of  stream 
function  on  pressure  surfaces  and  the 
nondivergent  wind  components  on  sigma  surfaces. 
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Fig.  18.  Flow  chart  diagramming  the  calculation  of  omega 
on  pressure  surfaces  and  the  divergent  wind 
components  on  sigma  surfaces. 
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Fig.  19.  Flow  chart  depicting  the  calculation  of  a. 
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calculated  on  sigma  surfaces.  Fig.  18  illustrates  how  the  divergent 
wind  components  are  obtained  on  sigma  surfaces.  In  Fig.  19,  we  show 
how  a  is  calculated.  Finally,  in  Fig.  20,  we  graphically  compare 
the  different  initialization  procedures.  Although  we  are  interested 
in  all  aspects  of  the  results  of  the  experiments,  we  will  be  primarily 
concerned  with  the  initial  precipitation  amounts. 

7.1  Initial  conditions 

In  the  three  forecast  experiments  that  were  balanced  (Forecasts  3, 
4,  and  5),  a  superadiabatic  lapse  rate  occurred  over  much  of  the 
surface  low  pressure  system  (the  eastern  Rockies  and  the  Plains 
states) .  In  each  initialization,  geopotential  is  defined  at  the 
surface  (sigma  level  7).  In  Forecast  3,  geopotential  is  calculated 
at  850  mb  while  in  Forecasts  4  and  5,  geopotential  is  calculated  at 
sigma  level  6^.  Temperatures  are  then  computed  for  the  bottom  layer 
from  the  hydrostatic  equation.  In  each  case,  the  lowest  layer  is 
about  500  m  thick.  One  explanation  for  the  superadiabatic  lapse 
rates  is  the  procedure  used  for  the  diagnosis  of  surface  pressure. 

The  initialization  procedure  uses  temperature  to  diagnose  surface 
pressure  from  sea  level  pressure  in  an  iterative  process  using  the 
hydrostatic  equation.  The  sea  level  pressure  observation  was 
originally  obtained  from  observed  station  pressure  using  the  standard 
atmospheric  lapse  rate.  Consequently,  if  the  actual  lapse  rate  is 
warmer  (or  colder)  than  the  standard  atmospheric  lapse  rate,  the 
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Fig.  20.  Flow  chart  illustrating  the  differences  in 
the  initializations  for  the  five  forecast 
experiments . 
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diagnosed  surface  pressure  will  be  higher  (or  lower)  than  the  observed 
surface  pressure.  In  this  case,  the  surface  pressure  was  under¬ 
estimated  in  the  vicinity  of  the  low  since  the  atmosphere  was  colder 
than  the  standard  atmospheric  lapse  rate  would  predict  for  the  850  mb 
to  surface  layer. 

Now  we  can  examine  typical  errors  in  the  temperature  at  the 
bottom  level  that  could  result  if  pg  or  the  geopotential  at  the  lowest 
calculation  level  are  in  error.  The  finite-difference  hydrostatic 
equation  for  p  =  0  can  be  written: 


In 


RT 


(7.1) 


where  T  applies  at  the  logarithmic  mean  pressure  of  p7  and  p^,  and 
the  subscripts  6%  and  7  refer  to  sigma  levels  64  and  7,  respectively. 
For  p^  =  1000  mb,  a  10  m  change  in  the  geopotential  at  level  64 
results  in  a  two-degree  change  in  the  mean  temperature  of  the  layer. 

A  3-mb  error  in  p^  (pg)  results  in  a  one-degree  change  in  the  mean 
temperature.  Keyser  (1978)  reported  that  for  a  925-mb  height  error 
of  10  m,  the  corresponding  error  in  temperature  for  the  925-mb  to 
surface  layer  is  3.8°C  for  SLP  *  1013  mb  and  4.4°C  for  SLP  »  1000  mb. 
Keyser  (personal  communication)  reported  the  occurrence  of  super- 
adiabatic  lapse  rates  (and  not  just  in  the  lowest  layer)  in  cases 
reported  by  Keyser  (1978)  and  Anthes  (1978).  Barker*  (personal 
communication)  reported  that  areas  with  a  superadiabatic  lapse  rate 
in  the  lowest  layer  occurred  after  balancing  the  initial  conditions 
* 

E.  H.  Barker,  1978,  Naval  Environmental  Prediction  Research  Facility, 

Monterey,  California. 
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for  a  global  synoptic-scale  primitive  equation  model.  Hence,  for  thin 
layers  especially,  the  hydrostatic  equation  is  very  sensitive  to  small 
errors  in  geopotential  or  surface  pressure. 

The  atmosphere  would  not  be  expected  to  have  a  super adiabatic 
lapse  rate  over  any  large  area.  Hence  the  forecast  model  should  not 
be  provided  with  these  erroneous  lapse  rates.  Keyser  and  Barker 
(personal  communication)  recommended  using  a  convective  adjustment 
procedure  at  the  end  of  the  initialization.  Although  the  balancing 
would  be  somewhat  altered,  the  result  is  more  realistic.  Also,  the 
presence  of  superadiabatic  layers  would  artificially  increase  the 
initial  precipitation  rate.  The  convective  adjustment  procedure 
selected  for  use  here  is  analogous  to  the  one  described  by  AW.  Their 
scheme  conserved  the  vertical  integral  of  internal  and  potential 
energy. 

The  initial  mean  values  of  wind  and  temperature  are  given  in 
Table  6.  The  temperatures  for  Forecasts  3,  4,  and  5  apply  after  the 
convective  adjustment  was  performed.  The  temperature  changes  produced 
by  the  convective  adjustment  procedure  were  several  degrees  at  levels 
5%  and  64  and  one  degree  at  levels  1%  and  24. 

Barker  (personal  communication)  found  that  the  number  of 
occurrences  of  superadiabatic  lapse  rates  dropped  dramatically  when 
a  new  three-dimensional  analysis  procedure  was  applied  to  the  raw 
data  before  balancing.  The  three-dimensional  analysis  scheme  replaced 
a  more  conventional  two-dimensional  system.  The  analysis  used  on 
the  12191175  data  was  a  two-dimensional  analysis.  This  suggests  that 
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Table  6.  Mean  values  of  Che  initial  temperature  and 
wind  fields  for  Forecasts  1  through  5. 


Forecast 

Number 

1 

2 

3 

4 

5 

Temperature 
Sigma  level 

00 

l*s 

234.1 

234.1 

235.0 

235.0 

235.0 

24 

254.1 

254.1 

252.6 

255.1 

255.0 

34 

262.9 

262.9 

261.3 

263.4 

263.3 

44 

268.7 

268.7 

268.8 

268.1 

268.0 

54 

274.1 

274.1 

274.3 

274.3 

274.3 

64 

279.7 

279.7 

280.2 

277.5 

V 

277.7 

u  component 
Sigma  level 

(m  s  1) 

14 

23.1 

23.1 

23.1 

23.1 

23.1 

24 

15.8 

15.9 

15.9 

15.9 

15.7 

34 

12.2 

12.2 

12.2 

12.2 

12.1 

44 

9.82 

9.84 

9.84 

9.84 

9.75 

54 

7.65 

7.67 

7.67 

7.67 

7.54 

64 

4.10 

4.12 

4.12 

4.12 

4.10 

v  component 
Sigma  level 

(m  s  1) 

14 

8.20 

8.21 

8.21 

8.21 

8.44 

24 

7.20 

7.18 

7.18 

7.18 

7.27 

34 

7.17 

7.17 

7.17 

7.17 

7.18 

44 

7.58 

7.59 

7.59 

7.59 

7.58 

54 

8.20 

8.19 

8.19 

8.19 

8.19 

64 

8.66 

8.65 

8.65 

8.65 

8.72 

an  analysis  procedure  with  some  vertical  consistency  constraint  should 
be  used  in  the  initialization  of  the  PSU  model  to  help  eliminate 
geopotential  errors. 

The  same  surface  pressure  field  was  used  in  all  five  experiments. 
Because  the  temperature  field  at  the  lowest  level  is  different  for  the 
unbalanced  experiments  and  each  balanced  experiment,  the  initial  sea 
level  pressure  fields  are  slightly  different.  Fig.  21a  is  the  sea 
level  pressure  field  for  Forecast  5.  Fig.  21b  is  the  balanced  500-mb 
temperature  field  and  Fig.  21c  is  the  500-mb  height  field  for 
Forecast  5.  The  lowest  height  in  the  trough  over  the  Rocky  Mountains 
is  15  m  less  than  in  Forecast  3.  Fig.  21d  is  the  initial  500-mb  wind 
field  for  Forecast  5.  The  nondivergent  wind  fields  for  Forecasts  2, 

3,  and  4  are  identical. 

One  method  of  comparison  of  the  initializations  involves  the 
computation  of  RMS  temperature  difference  between  them.  The  RMS 
temperature  differences  between  the  initial  conditions  of  the 
experiments  are  presented  in  Table  7.  The  values  are  not  unusual 
except  for  the  level  6*5  differences  between  the  balanced  experiments 
and  the  unbalanced  experiments.  An  examination  of  the  temperature 
difference  field  between  Forecasts  1  and  5  reveals  temperature 
differences  of  up  to  10°C  in  the  ridge  over  the  southeastern  United 
States  at  level  64*  In  Forecasts  4  and  5,  geopotentials  were 
calculated  on  sigma  surfaces  and  the  height  of  the  level  6*5  sigma 
surface  is  as  much  as  35  m  lower  than  the  same  sigma  surface  for 
Forecasts  1  and  3.  This  height  difference  is  due  to  the  mesoscale 
sigma-surface  balance  equation  used  for  Forecasts  4  and  5.  The  35-m 


Fig.  21.  Initial  fields  for  Forecast 


21.  (Continued). 


130 


height  difference  changes  the  depth  of  the  layer  between  sigma  levels 
6%  and  7  from  about  470  m  to  about  440  m.  This  represents  a  6  percent 
change  in  thickness  and  corresponds  to  a  6  percent  reduction  in  the 
mean  temperature  of  the  lowest  layer  of  about  20°C.  This  mean 
temperature  applies  at  about  sigma  level  6  3/4.  When  the  mean 
temperature  is  used  to  obtain  the  temperature  at  sigma  level  64,  the 
resultant  level  64  temperature  is  up  to  10°C  lower  than  the  unbalanced 
forecasts.  We  will  see  that  this  erroneous  temperature  difference 
adversely  affects  the  RMS  temperature  errors  at  the  lower  levels. 

The  geopotential  errors  in  the  vicinity  of  the  high  are  an  additional 
indication  that  an  analysis  procedure  with  some  vertical  constraint 
should  be  used  in  the  initialization  of  the  PSU  model.  That  is, 
for  a  thin  layer,  a  relatively  small  change  in  geopotential  can  make 
a  relatively  large  change  in  the  mean  temperature  of  the  layer. 

7.2  The  model  forecasts 

In  this  section  the  forecast  results  will  be  compared  and 
contrasted  with  the  exception  of  the  precipitation  predictions,  which 
will  be  discussed  in  Section  7.3.  For  each  forecast,  noise 
characteristics  and  the  12-hour  forecasts  of  sea  level  pressure, 
temperature,  and  wind  will  be  discussed. 
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7.2.1  Gravity-inertia  wave  characteristics  of  the  forecasts 


Gravity- inertia  waves  in  a  forecast  are  sometimes  referred  to 
as  "noise".  As  pointed  out  in  Chapter  1,  these  waves  can  destroy 
the  meteorological  portion  of  a  forecast  if  care  is  not  exercised 
during  the  initialization  and  the  forecast  itself.  However,  gravity 
waves,  particularly  internal  gravity  waves,  are  responsible  for  the 
geostrophic  adjustment  process.  Noise  is  perhaps  an  unfortunate 
choice  of  words  as  it  implies  that  all  gravity  waves  or  vertical 
motions  are  undesirable.  The  undesirable  vertical  motions  arise 
principally  from  the  time-dependent  boundary  conditions  and  aliasing 
and  are  generally  of  smaller  scale  than  the  meteorologically  significant 
vertical  motions.  The  undesirable  vertical  motions  are  minimized  with 
the  enhanced  horizontal  diffusion  scheme. 


Most  investigators  have  previously  used  a  time  series  of  the 


vertically  integrated  mass  divergence  as  an  indicator  of  noise.  That 

^P*  3— p^  ^P* 

is,  |-r— |  or  | — =—  j  have  been  used.  Sundqvist  (1975)  plotted  |tt— | 

3t  g  c  3t 

summed  over  the  domain  versus  time.  His  curves  showed  an  increase  for 


the  first  two  hours  after  initialization  and  then  displayed  a  gradual 
decline.  —9 — 


I  *  I 

Bleck  (1977)  preferred  to  use  | — =— |  .  He  stated  that  gravity 

3t  3D 

waves  appeared  most  clearly  in  the  —  term  of  the  divergence 

m  O  t 

1  ^nP* 

equation  and,  therefore,  | — —  j  should  be  used  to  indicate  external 


gravity-wave  noise.  Since 


_*|  is  related  to  —  dp  by  the 
?  — X  _  3t 


,3t2  ,=2p. 


I  0  r  *  I  S 

continuity  equation,  Bleck  used _  j - 1  as  a  measure  of  general 

3ZP*  3t2 

noisiness.  Bleck's  plots  of  j - 1  showed  a  sharp  decline  in  general 

3t 
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noisiness  for  the  first  three  hours  followed  by  a  gradual  decline.  His 


initial  conditions  were  unbalanced  so  we  should  expect  a  high  initial 

d2p  * 

noise  level.  Bleck’s  model  with  terrain  produced  a  value  for  -  of 

a*2 

-7  -2  at 

about  5  x  10  mb  s _ at  the  end  of  a  12-hour  forecast  on  an  85-km  mesh. 

3p*  32p^ 

WAM  reported  |-^£— |  and  | - 1  for  the  PSU  model  applied  to  another 

3t 

data  set.  A  time  smoother  has  since  been  incorporated  into  the  PSU 


model  and  the  noise  characteristics  should  now  be  somewhat  different 
than  those  in  WAM.  Anthes  (1978)  used  a  version  of  the  model  employed 

here,  on  a  60-km  grid  mesh  and  another  data  set.  He  reported  mean 

-4  -1  -7  -2  ,3P*. 

values  of  about  5  x  10  mb  s  and  5  x  10  mb  s  for  |-r—  |  and 
0  dt 

i3  p*i 

I - 1 ,  respectively,  at  the  end  of  the  first  12  h  of  a  24-hour 

3t2 

forecast.  _  — r — 

3p*  3  P* 

Figs.  22  and  23  are  the  plots  of  |-r—  |  and  | — =— |  ,  respectively, 

3t  3t 

versus  time  for  Forecasts  1  through  5.  In  Fig.  22,  the  curves  increase 

slightly  for  about  an  hour  and  then  decrease  slowly  through  the  rest 

3P* 

of  the  forecast  period.  Forecasts  1  and  2  have  similar  |— -|  plots 

d  t 

and  are  not  significantly  different  in  this  respect.  Likewise, 
Forecasts  3,  4,  and  5  have  similar  graphs.  The  forecasts  started 


from  balanced  initial  conditions  (Forecasts  3,  4,  and  5)  do  indicate 

3P* 

.  I  75  1  _ 


a  lower  noise  level  as  measured  by 


than  the  unbalanced 


experiments,  1  and  2.  In  Fig.  23,  the  plots  decrease  rapidly  for 
three  hours  and  then  decrease  more  slowly  for  the  rest  of  the  forecast. 
The  forecasts  with  unbalanced  initial  conditions  again  possess  a  higher 
noise  level  than  the  balanced  experiments.  Forecast  1  has  a 
significantly  higher  noise  level.  As  in  Fig.  22,  Forecasts  3,  4,  and 
5  are  not  substantially  different.  For  Forecast  5,  this  implies  that 


l  /V'  *  v** v~’  -*  - 


1  through  5. 
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the  added  divergence  does  not  generate  significant  additional  noise 
in  the  forecast.  That  is,  very  little  of  the  divergence  is  partitioned 
into  gravity-inertia  wave  noise.  In  both  Figs.  22  and  23,  we  detect 
no  difference  in  noise  level  for  the  balanced  experiments. 

The  results  in  Figs.  22  and  23  are  similar  to  those  of  Bleck 

i s2P* i 

(1977)  and  Anthes  (1978).  Bleck' s  value  of  | - 1  for  unbalanced 

3t2 

initial  conditions  at  12  h  was  the  same  as  that  in  Fig.  23  for  Forecast 


3P* 

1.  Anthes  reported  slightly  higher  values  of  I'^H  and 


32p*| 

V 


at  the 


end  of  12  h  than  for  the  balanced  experiments  in  Figs.  22  and  23, 

respectively.  This  was  probably  due  to  his  use  of  a  60-km  grid  mesh, 

compared  to  the  120-km  mesh  used  here.  In  another  forecast  experiment 

not  reported  on  here,  the  parameter  which  controls  the  strength  of 

the  enhanced  horizontal  diffusion  was  increased  by  a  factor  of  five. 

Although  the  precipitation  amounts  on  the  domain  interior  were  not 

3P* 

affected,  the  noise  level  as  measured  by  |t— — |  was  reduced  by  about  a 


factor  of  three.  Hence  the  temporal  variation  of  surface  pressure 
depends  to  some  degree  on  the  horizontal  diffusion  employed  near  the 
boundaries.  This  fact  alone  could  account  for  the  slight  differences 
between  the  noise  statistics  reported  by  Bleck  and  Anthes  and  those 
reported  here. 

Sundqvist  (1975)  reported  a  significant  reduction  of  noise  level 
when  the  mass  and  momentum  variables  were  balanced  on  sigma  surfaces 
for  a  synoptic-scale  model.  He  also  reported  a  higher  noise  level 
for  sigma-surface  balancing  when  mountains  were  used  as  the  lower 
boundary.  The  terrain  field  used  by  Sundqvist  was  smoothed  such  that 
4Ax  and  shorter  wavelengths  were  removed.  No  smoothing  was  applied 
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to  the  terrain  field  used  in  these  experiments.  Figs.  22  and  23  show 
no  noise  level  reduction  for  balancing  on  sigma  surfaces.  Perhaps 
the  mesoscale  grid  mesh,  limited  domain,  and  terrain  field  employed 
did  not  permit  a  significant  noise  reduction  for  Forecast  4,  as 
compared  to  Forecast  3. 

0kland  (1970)  studied  a  relatively  simple  two-level  baroclinic 

aP* 

model.  He  plotted  RMS  — —  versus  time  for  nondivergent  and  divergent 

initial  conditions  (divergent  initial  conditions  consisted  simply  of 

the  model  pred^  tion  after  12  hours  of  integration;  that  is,  the 

forecast  was  started  12  hours  previously  from  different  data) . 

3p* 

0kland's  RMS  — —  curve  was  noisier  for  nondivergent  initial  conditions. 

3P* 

As  he  pointed  out,  — —  was  primarily  a  measure  of  external  gravity  wave 
a  t 

activity.  Therefore,  the  forecast  started  from  nondivergent  initial 
conditions  had  more  external  gravity  wave  activity.  0kland  also 
plotted  RMS  omega  values  for  a  level  near  500  mb.  He  pointed  out 
that,  because  internal  gravity  modes  have  higher  vertical  velocities 
associated  with  them  than  external  modes,  the  RMS  omega  curve  will 
reflect  primarily  the  internal  gravity  waves.  The  RMS  omega  curve 
given  by  0kland  for  the  divergent  initial  conditions  was  relatively 
constant  indicating  that  the  internal  modes  were  fully  developed  and 
that  the  model  variables  were  approximately  in  balance.  In  0kland's 
nondivergent  case,  the  RMS  omega  curve  began  at  zero,  reached  the 
divergent  initial  condition  omega  curve  in  3  to  4  hours,  continued  to 
rise  to  a  maximum  in  6  to  7  hours,  and  then  decreased  to  a  minimum 


at  11  to  13  hours. 
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Here  precipitation  rates  are  the  primary  interest.  External 
gravity  waves  usually  have  a  relatively  small  effect  on  precipitation. 
Precipitation  arises  mostly  from  the  slower  internal  gravity  modes. 
Therefore,  the  500-mb  RMS  omega  graphs  versus  time  will  be  presented 
for  Forecasts  3,  4,  and  5.  Although  the  previous  measures  of  gravity- 
wave  intensity  (Figs.  22  and  23)  were  very  similar  for  the  three 
forecasts,  the  precipitation  amounts,  as  shall  be  shown  later,  are 
significantly  different.  The  RMS  omega  graphs  will  be  related  to  the 
forecast  precipitation. 

Fig.  24  contains  the  500-mb  RMS  omega  plots  for  Forecasts  3,  4, 
and  5  (the  plots  for  1  and  2  are  very  similar  to  3  and  4) .  The  curve 
for  5  begins  at  a  higher  RMS  omega  value  and,  after  5  hours,  reaches 
a  maximum  value  lower  than  that  of  the  nondivergent  experiments.  A 
strong  similarity  exists  between  Fig.  24  and  0kland's  graph  of  the 
RMS  omega.  Forecasts  2,  3,  and  4  begin  with  no  divergence.  As  the 
geostrophic  adjustment  process  proceeds  and  those  experiments  begin 
to  develop  a  divergent  component,  they  actually  "overshoot"  the 
quasi-equilibrium  value  and  develop  a  larger  divergent  component  than 
Forecast  5  from  4  to  6  hours  into  the  forecast.  The  500-mb  RMS 
omega  curve  for  Forecast  5  indicates  the  rate  of  precipitation  will 
rise  slowly  during  the  first  hours  of  the  forecast  and  then  remain 
relatively  constant. 

Forecasts  1  and  2  possess  higher  noise  levels  than  Forecasts  3, 

32p*i 

4,  and  5.  Forecast  1  is  particularly  noisy  as  measured  by  ] - ]  . 

3 1 

The  noise  characteristic;-  of  Forecasts  3,  4,  and  5  as  measured  by 
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graphs  of  *r—  and  | — r— I  are  not  significantly  different.  That  is, 
3t  3C2 

no  noise  reduction  resulted  by  balancing  on  sigma  surfaces.  However, 
the  500-mb  RMS  omega  curves  are  significantly  different  between 


Forecast  5  and  the  nondivergent  forecast  experiments. 


7.2.2  The  forecasts  of  sea  level  pressure,  temperature,  and  wind 


Table  8  contains  the  RMS  errors  in  the  forecasts  for  sea  level 
pressure  and  for  temperature  and  the  wind  components  at  400,  500,  and 
700  mb.  No  significant  difference  in  the  forecasts  can  be  detected 
from  these  statistics  although  it  is  interesting  to  note  that  the 
unbalanced  forecasts  (Forecasts  1  and  2)  had  the  lowest  RMS  temperature 
errors.  Table  8  also  contains  the  SI  scores  for  sea  level  pressure 
for  each  forecast.  The  SI  score  was  developed  by  Teweles  and  Wobus 
(1954).  It  is  an  objective  measure  of  the  forecast  skill  of  the  sea 
level  pressure  forecast  that  is  simple  to  compute  and  has  been 
reported  often  for  other  models.  The  SI  score  for  sea  level  pressure 
relates  horizontal  differences  in  the  forecast  sea  level  pressure  to 
the  observed  differences.  The  SI  score  is  computed  from 


SI  -  100 


Z  AP,  -  AP 
,  ,  1  f  o 

1-1 _ 

12 

Z  G, 


where  AP,  and  AP  are  the  differences  in  forecast  and  observed  sea 
f  o 

level  pressure,  respectively,  and  is  the  maximum  of  AP^  and  APq 
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for  a  given  pair  of  points.  The  differences  in  sea  level  pressure 
were  calculated  for  all  adjacent  points  in  the  set  of  nine  points 
listed  in  Table  9,  which  gives  a  total  of  12  combinations  of  points. 
The  "observed"  sea  level  pressure  differences  were  calculated  from 
the  analysis  of  sea  level  pressure  given  in  Fig.  10a.  An  SI  score 
of  20  is  considered  perfect  for  practical  purposes  while  a  score  of 
70  is  nearly  worthless. 

With  a  slightly  different  version  of  the  model  used  here,  Anthes 
(1978)  reported  an  average  SI  score  at  12  hours  of  33.8  for  32  cases 
over  the  US  and  Europe.  The  FNWC  hemispheric  model  achieved  an  SI 
score  at  12  hours  of  37.2.  These  scores  demonstrated  considerable 
skill  at  12  hours. 

The  SI  scores  reported  here  compare  well  with  those  of  Anthes 
and  those  at  FNWC.  It  is  interesting  to  note  that  the  best  SI  scores 
here  are  for  the  forecasts  from  unbalanced  initial  conditions 
(Forecasts  1  and  2).  The  SI  score  of  30.2  for  the  forecast  from 
observed  initial  conditions  was  significantly  better  than  the  other 
experiments.  One  conclusion  from  Table  8  is  that  divergent 
initialization  did  not  improve  the  SI  score  for  sea  level  pressure 
or  the  RMS  errors  in  the  temperature  and  wind  fields. 

0kland  (1970)  obtained  an  analytic  solution  for  a  simple 
baroclinic  model.  From  this  solution,  he  concluded  that  the  final 
adjusted  state  is  independent  of  the  initial  divergence.  Therefore, 
the  RMS  errors  for  the  experiments  in  Table  f  should  not  be  expected 
to  be  much  different.  Divergent  initialization  will,  however, 
influence  the  initial  precipitation  rate. 
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!  9. 

Indices  of 

the  points 

used  in 

computing  SI 

scores. 

i  denotes 

south-north 

direction;  j  denotes 

west-east 

direction, 
is  north 

Point  (1, 
latitude ; 

1)  is  the 
X  is  west 

lower  left  i 
longitude. 

corner. 

Point 

i 

i 

i 

\ 

1 

5 

10 

29.4 

104.2 

2 

5 

20 

29.5 

91.8 

3 

5 

30 

28.1 

79.7 

4 

15 

10 

40.3 

106.0 

5 

15 

20 

40.8 

91.3 

6 

15 

30 

39.0 

76.7 

7 

25 

10 

51.5 

108.8 

8 

25 

20 

52.1 

90.6 

9 

25 

30 

50.2 

72.5 
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The  sea  level  pressure  forecasts  for  the  experiments  are  similar. 
Fig.  25a  is  the  sea  level  pressure  forecast  for  Forecast  5.  The  main 
difference  between  Fig.  25a  and  the  nondivergent  experiments  is  that 
the  low  in  the  nondivergent  forecasts  is  1009  mb.  The  observed 
central  pressure  of  the  low  was  1004  mb  (Fig.  10a). 

The  500-mb  temperature  fields  for  the  forecasts  are  again  similar. 
Fig.  25d  is  the  500-mb  temperature  field  for  Forecast  5.  None  of 
the  forecasts  extended  the  cold  air  over  the  Rocky  Mountains  far 
enough  to  the  south  (compare  with  Fig.  lOd) .  The  intensity  of  the 
trough  is  underforecast. 

Fig.  25b  is  the  500-mb  height  forecast  for  Forecast  5.  The 
contour  pattern  is  similar  to  the  other  forecasts.  The  lowest  height 
over  Colorado  for  Forecast  5  was  12  m  lower  than  in  Forecast  3.  The 
observed  height  was  5472  m  over  the  Colorado-Kansas  border  (see 
Fig.  10b). 

Fig.  25c  is  the  500-mb  wind  forecast  for  Forecast  5.  The 
forecasts  are  similar.  The  forecast  agrees  reasonably  well  with  the 
observed  wind  field  (Fig.  10c)  with  the  exception  of  the  wind  maximum 
near  El  Paso.  The  forecast  only  slightly  increased  the  core  speed 
of  the  jet  and  moved  the  core  northward  a  small  distance.  However, 
the  core  speed  of  the  jet  streak  actually  increased  by  15.4  m  s  ^ 
and  the  core  moved  to  a  point  near  Dallas.  As  with  the  temperature 
and  height  fields,  the  wind  forecast  reflects  the  fact  that  the  model 
did  not  forecast  the  rapid  intensification  of  the  trough  or  enough 
deepening  of  the  surface  low. 
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7.3  The  precipitation  forecasts 

In  this  section  the  precipitation  forecasts  will  be  compared. 

The  precipitation  forecasts  will  be  viewed  in  several  ways.  First  the 
forecasts  will  be  "scored".  Then  the  precipitation  amounts  themselves 
will  be  examined. 

7.3.1  Scoring  the  precipitation  forecasts 

Anthes  (unpublished)  developed  an  objective  procedure  for  the 

"scoring"  of  precipitation  fields.  If  the  forecast  precipitation  field 

is  one  grid  increment  displaced  from  the  observed  precipitation  field, 

then  the  forecast  still  contains  useful  information.  The  objective 

procedure  that  will  be  used,  scores  a  precipitation  forecast  in  a 

quantitative  way  to  determine  if  a  shift  in  the  forecast  field  with 

respect  to  the  observed  field  improves  the  forecast  when  compared  with 

the  observations.  The  complete  objective  procedure  is  presented  in 

Appendix  3.  Therefore,  it  is  discussed  only  briefly  here.  A 

correlation  coefficient  is  calculated  for  each  shift  in  position  of  the 

forecast  field  with  respect  to  the  observed  field.  An  example  matrix 

of  correlation  coefficients  c^  ^  for  a  shift  in  each  direction  of  two 

grid  increments  is  given  in  Appendix  3.  For  example,  Cq  q  is  the 

correlation  coefficient  for  no  shift  and  c..  _  is  the  correlation 

1 ,  u 

coefficient  for  a  shift  one  grid  increment  to  the  right.  Good  forecasts 
would  have  a  large  maximum  c^  ^  (1.0  is  perfect)  and  low  values  for  l 
and  k  for  the  maximum  correlation  coefficient  (0  and  0  are  perfect). 


Figs.  26a,  26b,  and  26c  contain  the  correlation  coefficient 
matrices  for  the  first  three  hours  of  the  forecast  for  Forecasts  3, 

4,  and  5,  respectively  (Forecasts  1  and  2  are  similar  to  Forecast  3). 
Note  that  for  the  nondivergent  forecasts  (3  and  4),  the  highest 
correlation  coefficient  is  displaced  one  grid  increment  from  the 
center.  For  Forecast  5,  the  highest  correlation  coefficient  occurs 
with  no  spatial  offset.  Also  the  maximum  correlation  coefficient 
for  Forecast  5  (c^  ^  =  0.75)  is  higher  than  the  nondivergent  forecasts 
(Cq  ^  =  0.55  for  Forecast  3  and  cg  ^  =  0.66  for  Forecast  4).  That  is, 
the  forecast  experiment  started  with  a  divergent  initialization 
produced  the  best  scores. 

Figs.  27a,  27b,  and  27c  contain  the  correlation  coefficient 
matrices  for  the  entire  12-hour  period  for  Forecasts  3,  4,  and  5. 

The  maximum  correlation  coefficients  for  the  three  forecasts  were 
0.77,  0.75,  and  0.76  for  Forecasts  3,  4,  and  5,  respectively.  For 
the  entire  forecast  period,  all  three  forecast  experiments  had  the 
highest  correlation  coefficient  with  no  spatial  offset  of  the  forecast 
and  observed  precipitation  fields.  Therefore,  the  objective  procedure 
detects  no  differences  in  the  12-hour  precipitation  forecasts  when 
divergence  is  included  in  the  initial  conditions. 

7.3.2  The  total  precipitation  amounts  and  the  precipitation  rates 

Fig.  28  contains  the  total  precipitation  amounts  for  each  hour 
through  the  forecast  period  for  each  experiment  as  well  as  the  total 
depth  of  observed  precipitation.  The  amounts  were  summed  for  the 
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Fig.  26.  Correlation  coefficient  matrices  from  the 
objective  precipitation  scoring  procedure 
for  the  first  three  hours  of  the  forecast. 
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•  27.  Correlation  coefficient  matrices  from  the 
objective  precipitation  scoring  procedure 
for  the  entire  12-hour  forecast  period. 
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entire  domain  in  order  to  compare  the  total  depth  of  precipitation 
for  each  forecast  with  the  observed  precipitation.  Multiplying  each 
depth  by  the  area  of  a  grid  square  would  yield  precipitation  volume. 

The  total  amount  of  precipitation  forecast  by  the  model  is  a  useful 
statistic  because  the  forecasts  over  the  whole  domain  can  be  compared. 
Similarly,  Figs.  29a  and  29b  present  the  total  precipitation  depth 
forecast  by  the  convective  and  nonconvective  precipitation  parameteriza- 
tions,  respectively. 

For  the  nondivergent  experiments,  a  slight  decrease  in  the  total 
precipitation  occurs  from  hour  1  to  hour  2  followed  by  an  increase 
each  hour  to  the  maximum  amounts  in  hours  7,  8,  and  9.  Then  the 
amounts  generally  decrease  for  each  hour  through  the  end  of  the 
forecast.  The  behavior  of  the  convective  and  nonconvective  precipitation 
follows  the  same  general  pattern.  However,  Forecasts  1  and  2  have 
consistently  smaller  nonconvective  precipitation  amounts  after  hour  3. 

In  contrast,  the  total  precipitation  for  the  forecast  from  a  divergent 
initialization  (Forecast  5)  is  more  uniform  throughout  the  forecast 
period.  The  precipitation  from  Forecast  5  does  reach  maximum  values 
in  hours  7,  8,  and  9,  and  these  maximum  values  are  approximately 
equal  to  Forecasts  3  and  4  for  those  hours.  However,  the  maximum 
values  are  only  about  20  percent  greater  than  the  first  and  last 
forecast  hours.  The  observed  precipitation  is  about  uniform  throughout 
the  forecast  period  although  there  is  an  increase  in  the  last  three 
hours.  Recall  the  500-mb  RMS  omega  curves  of  Fig.  24.  The  relatively 
constant  curve  for  Forecast  5  is  reflected  in  the  relatively  constant 
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Fig.  29.  Total  precipitation  depth  (cm)  for  the  entire 
domain  for  the  convective  and  nonconvective 
precipitation  parameterizations  for  Forecasts 
1  through  5. 
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precipitation  amounts  for  each  hour.  The  greater  variability  of  the 
curves  for  Forecasts  3  and  4  is  reflected  in  the  greater  variability 
in  the  precipitation  amounts  by  hour  for  those  experiments.  Note 
that  the  maxima  in  the  RMS  omega  curves  for  Forecasts  3  and  4  occur 
at  hour  5  while  the  maxima  in  their  precipitation  amounts  occur  later 
in  the  forecast  period.  Therefore,  the  500-mb  RMS  omega  curves  are 
not  perfectly  correlated  with  the  precipitation  amounts.  Instead, 
it  is  suggested  that  a  greater  variability  in  the  RMS  omega  curves 
implies  a  greater  variability  in  the  rainfall  amounts,  but  that  the 
maxima  in  the  two  do  not  necessarily  coincide. 

The  12-hour  precipitation  depths  summed  for  every  grid  point 
over  the  domain  for  the  experiments  are  55.2  cm,  50.5  cm,  58.8  cm, 

65.6  cm,  and  72.9  cm  for  Forecasts  1  through  5,  respectively.  Note 
that,  for  the  entire  period,  Forecast  5  produced  32,  44,  24,  and  11 
percent  more  precipitation  than  Forecasts  1,  2,  3,  and  4,  respectively. 
However,  the  12-hour  precipitation  total  for  Forecast  5  is  only  58.1 
percent  of  the  observed  precipitation  total.  The  12-hour  observed 
precipitation  is  really  an  approximation  because  precipitation  is 
sampled  only  at  a  few  points.  Also,  there  is  evidence  of  convective 
activity  in  Fig.  lib  (small-scale  peaks  in  precipitation  amounts). 

A  model  with  a  grid  increment  of  120  km  cannot  be  expected  to  forecast 
small  (e.g.,  30  km  by  30  km)  areas  of  convective  activity.  For  this 
reason,  the  model  should  not  be  expected  to  forecast  100  percent  of 
the  observed  point  values  of  precipitation,  when  these  values  are 
not  representative  of  grid-square  averages.  In  spite  of  this,  the 
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precipitation  forecast  for  the  divergent  initialization  is  a 
significant  improvement  over  the  other  forecasts. 

For  all  the  experiments,  the  total  precipitation  for  hours  7 
through  12  is  about  the  same.  Therefore,  the  difference  in  the 
12-hour  forecast  amounts  is  due  almost  entirely  to  the  first  6  hours 
of  the  forecasts.  Additionally,  most  of  the  difference  in  the  first 
6  hours  occurs  in  the  first  3  of  those  6  hours.  In  the  first  3  hours. 
Forecast  5  produced  89,  118,  83,  and  56  percent  more  precipitation 
than  Forecasts  1,  2,  3,  and  4,  respectively.  This  precipitation 
increase  is  directly  attributable  to  the  divergent  initialization. 

When  nondivergent  initial  conditions  are  used,  the  model  did 
produce  about  the  same  amount  of  precipitation  in  the  last  half  of 
the  forecast  as  it  did  with  divergent  initial  conditions.  The  non¬ 
divergent  forecasts  did  not  "overshoot".  That  is,  they  did  not 
precipitate  more  than  the  divergent  forecast  later  in  the  12-hour 
period.  If  the  nondivergent  forecasts  did  overshoot,  they  would  begin 
to  catch  up  to  the  total  12-hour  amount  produced  by  the  divergent 
initialization.  Therefore,  the  precipitation  lost  at  the  beginning 
of  a  forecast  started  from  a  nondivergent  initialization  is  not 
recovered . 

For  Forecast  5,  the  initial  slight  increase  in  the  500-mb  RMS 
omega  curve  in  Fig.  24  and  the  more  realistic  precipitation  amounts 
in  the  first  4  hours  of  the  forecast  are  the  best  indicators  that  the 
model  did  indeed  remember  the  initial  divergence. 
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The  usefulness  of  divergent  initialization  is  not  just  the 
additional  precipitation  in  the  first  6  hours  of  the  forecast  but  also 
the  relatively  constant  precipitation  rate.  Although  precipitation 
does  not  necessarily  occur  at  an  almost  constant  rate,  it  did  not 
occur  in  the  pattern  of  amounts  shown  by  the  nondivergent  forecasts 
either.  That  is,  the  amounts  in  Fig.  28  for  the  nondivergent  experiments 
reflect  the  anomolous  divergence  associated  with  the  geostrophic 
adjustment  process  in  the  model.  The  divergent  component  gradually 
builds  up,  overshoots  the  balanced  state,  and  then  gradually  oscillates 
about  the  balanced  state.  This  is  the  principle  cause  of  the 
oscillation  in  the  precipitation  amounts  for  the  nondivergent  forecasts. 
However,  in  the  divergent  initialization  experiment,  the  model  was 
closer  to  an  initial  balance  between  the  nondivergent  and  divergent 
components.  Hence  there  is  less  tendency  for  the  model  to  reflect 
temporal  maxima  and  minima  in  the  mean  precipitation  amounts  that  are 
as  extreme  as  for  the  nondivergent  forecasts. 

The  precipitation  fields  for  Forecasts  3,  4,  and  5  will  now  be 
compared  with  each  other  and  with  the  observed  precipitation. 

Figs.  30b,  30c,  and  30d  are  the  forecast  precipitation  amounts  for  the 
first  three  hours  for  Forecasts  3,  4,  and  5,  respectively.  Figs.  31b, 
31c,  and  31d  are  the  forecast  precipitation  amounts  for  the  12-hour 
forecasts  for  Forecasts  3,  4,  and  5,  respectively.  Figs.  30a  and 
31a  are  the  observed  precipitation  amounts  for  the  first  three  hours 
and  the  entire  12-hour  period,  respectively.  There  is  some  evidence 
of  noise  in  all  the  forecast  precipitation  figures.  For  example. 


31.  Observed  precipitation  for  the  entire  12-hour  forecast  period  (12191175  to  00201175)  and 
the  forecast  precipitation  for  the  same  period  for  Forecasts  3,  4,  and  5. 


(Continued) . 
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the  small  area  of  precipitation  in  the  northwest  corner  of  the  domain 
is  due  to  boundary-generated  noise  and  is  not  meteorological. 

In  comparing  Figs.  30b,  30c,  and  30d,  note  that  the  0.25-cm 
contour  in  Fig.  30d  covers  a  larger  area  than  in  Figs.  30b  and  30c. 

It  is  apparent  that  more  precipitation  occurred  in  the  first  three 
hours  of  the  forecast  with  divergent  initial  conditions.  For  the 
12-hour  forecast  precipitation  amounts,  the  0.25-cm  contour  in 
Figs.  31b,  31c,  and  31d  covers  about  the  same  area.  Again,  the 
highest  amounts  appear  inside  the  0.25-cm  contour  of  Fig.  31d. 

Since  the  model  did  produce  an  improved  precipitation  forecast 
with  the  divergent  initialization  used  in  Forecast  5,  especially  in 
the  first  few  hours,  the  initial  divergence  was  "remembered"  or 
retained  by  the  model.  That  is,  the  initial  divergence  was  not 
dissipated  by  internal  gravity  waves  before  it  could  be  supported 
by  precipitation  and  the  associated  release  of  latent  heat.  Since 
the  precipitation  occurred  in  the  correct  locations,  the  divergent 
wind  fields  and  consequently  the  diagnosed  omega  field  were 
sufficiently  accurate  to  be  useful  in  divergent  initialization. 

In  order  to  show  that  the  model  did  retain  the  initial  divergence, 
the  initial  precipitation  rates  at  several  locations  will  be 
examined.  These  locations  are  plotted  in  Fig.  32.  The  points  were 
selected  because  most  of  those  locations  had  significant  precipitation 
at  the  initial  time.  For  each  location  in  Fig.  32,  Fig.  33  contains 
the  forecast  precipitation  for  each  0.1-hour  interval  for  the  first 
1.6  hours  of  the  forecast  period.  Fig.  33  also  contains  the  "observed" 
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Fig.  32.  Locations  at  which  the  initial 
precipitation  was  examined. 
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Fig.  33.  Initial  precipitation  rates  (cm  h  at  the  locations 
in  Fig.  32  for  the  first  1.6  hours  of  the  forecast 
from  Forecast  5.  The  dashed  line  represents  the 
"observed"  initial  rain  rate  as  depicted  in  Fig.  13. 
The  solid  line  is  the  precipitation  amount  for  each 
two  time-step  period  (0.1  hour). 
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Fig.  33.  (Continued). 


initial  precipitation  rate.  The  observed  rate  was  taken  from  the 
same  four-hour  average  precipitation  rate  that  was  used  for  the 
diabatic  term  of  the  omega  equation  (Fig.  13).  The  purpose  of  Fig.  33 
is  to  show  how  well  the  forecast  initial  precipitation  rate  for 
Forecast  5  matches  the  observed  initial  precipitation  rate.  At  most 
of  the  points,  the  forecast  initial  precipitation  rate  is  near  the 
observed  rate.  Also,  the  forecast  rate  is  relatively  constant 
indicating  that  the  initial  divergence  did  not  dissipate  or  migrate 
away  from  the  initial  precipitation  area.  The  worst  result  was  at 
location  5.  It  is  near  an  area  of  relatively  light  observed  initial 
precipitation.  Also,  the  initial  moisture  analysis  indicates  that 
location  5  is  on  the  edge  of  the  saturated  area.  Therefore,  the 
initial  divergence  probably  did  not  initiate  precipitation  fast  enough 
for  the  divergence  to  be  maintained.  Recall  that  saturated  initial 
conditions  are  essential  for  divergent  initialization  in  areas  of 
initial  observed  precipitation  because  immediate  latent  heat  release 
is  required  to  maintain  the  upward  motion  associated  with  strong 
convergence.  Otherwise,  the  initial  divergence  would  be  dissipated 
by  internal  gravity  waves. 

Fig.  33  demonstrates  that  the  initial  divergence  is  retained 
reasonably  well.  However,  it  is  possible  that  the  initial  divergence 
produced  precipitation  at  an  unexpected  location.  For  example,  the 
precipitation  that  was  expected  at  location  5  could  have  occurred 
to  the  south  or  southwest  of  location  5.  This  possibility  will  be 
examined  in  Fig.  34.  Figs.  34b  through  34o  contain  the  initial 
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forecast  rain  rates  (0.05  cm  h  ^  contour  interval)  for  each  0.1-hour 
period  for  the  first  1.4  hours  of  the  divergent  initialization 
experiment  for  that  portion  of  the  domain  which  had  the  highest 
observed  initial  precipitation  rate.  For  comparison  purposes.  Fig. 

34a  contains  the  observed  initial  precipitation  rate  (0.05  cm  h  ^ 
contour  interval)  as  well  as  the  five  locations  from  Fig.  32  that 
fall  within  the  area  of  heaviest  precipitation.  First,  note  that  the 
zero  contour  is  forecast  well  in  the  northwest,  east,  southeast,  and 
southwest  portions  of  the  region  portrayed.  Second,  the  highest 
initial  precipitation  rates  occur  on  too  small  a  scale  to  be  retained. 
Sufficient  vertical  motion  to  sustain  them  was  probably  not  diagnosed 
by  the  omega  equation  on  that  scale  either.  The  largest  forecast 
rain  amounts  that  did  occur  are  in  Fig.  34b.  These  could  be  due  in 
part  to  the  relatively  high  initial  noise  level  since  no  attempt 
was  made  to  obtain  a  divergent  wind  component  that  was  exactly 
compatible  with  the  nondivergent  wind  component.  Third,  although 
the  forecast  0.1  cm  h  ^  contour  covers  a  smaller  area  than  the 
observed  one,  note  that  it  consistently  remains  about  the  same  size. 

It  is  almost  exactly  the  same  size  in  Figs.  34d  through  34h  and  again 
in  Figs.  341  through  34o.  There  are  no  rapid  changes  in  the 
forecast  pattern  that  would  be  associated  with  the  dissipation  of 
the  initial  divergence.  Fig.  34  shows  that  the  divergent  initialization 
experiment  did  not  forecast  the  total  observed  amount  because  the 
highest  observed  precipitation  rates  were  not  duplicated  and  because 
the  area  of  heaviest  forecast  precipitation  was  not  as  large  as  the 
corresponding  observed  area. 


Figs.  33  and  34  demonstrate  conclusively  that  the  initial 
precipitation  rates  were  steady  and  duplicated  the  observed  initial 
rates  reasonably  well.  The  divergent  initialization  was  useful  in 
improving  the  initial  precipitation  rates.  The  following  are  possible 
reasons  why  the  divergent  initialization  experiment  did  not  duplicate 
the  observed  initial  precipitation  rate  exactly: 

(1)  The  forecast  model's  cumulus  parameterization,  which 
was  responsible  for  most  of  the  initial  precipitation, 
might  not  be  producing  precipitation  at  the  correct 
rate  for  the  given  profiles  of  specific  humidity  and 
the  dynamic  variables. 

(2)  The  moisture  analysis  may  not  have  been  saturated 
over  a  large  enough  area  to  support  the  necessary 
latent  heating  rate  implied  by  the  initial  vertical 
motion  field. 

(3)  The  initial  divergent  wind  component  was  not  completely 
compatible  with  the  nondivergent  component. 

(4)  Some  of  the  initial  divergence  may  have  dissipated 
before  latent  heating  had  time  to  support  it. 

(5)  The  heaviest  observed  precipitation  rate  occurred  on 
a  scale  too  small  for  the  omega  equation  and/or  the 
forecast  model  to  represent  on  a  120-km  grid  mesh. 

For  this  case  divergent  initialization  significantly  improved  the 
short-range  forecasting  of  precipitation  with  a  mesoscale  numerical 


weather  prediction  model  started  from  a  static  initialization. 

Greater  improvements  would  be  likely  for  cases  with  higher  precipitation 
rates  at  the  initial  time. 

7.4  Chapter  summary 


In  this  chapter  five  12-hour  forecasts  were  compared  in  which 
the  forecast  model  was  initialized  in  different  ways.  The  five 
initializations  were  presented  in  Fig.  20. 


The  external  gravity-wave  noise  characteristics  of  the  forecasts 
.8P*.  - 32p* , 

as  measured  by  ]  -r— — |  and  | _ |  were  not  greatly  different  although 

o  t  2 

the  forecasts  from  unbalanced  initial  conditions  contained  more  noise. 


These  measures  of  noise  levels  were  comparable  to  those  values 
reported  by  Anthes  (1978)  and  Bleck  (1977).  Sundqvist  (1975)  reported 
that  a  lower  noise  level  resulted  from  an  initialization  on  sigma 
surfaces.  We  did  not  obtain  a  reduced  noise  level  in  the  experiments 
balanced  on  sigma  surfaces.  This  was  probably  due  to  the  smaller 
grid  increment,  unsmoothed  terrain,  and  limited  domain  used  here. 

The  500-tnb  RMS  omega  curve  for  the  divergent  initialization  (Forecast  5) 
was  significantly  different  from  others  as  it  exhibited  less 
variability  during  the  12-hour  forecast.  The  internal  gravity-wave 
activity  in  the  forecast  from  the  divergent  initialization  was 
therefore  more  uniform  and  the  model  was  initially  closer  to  a 
"balanced"  state.  The  500-mb  RMS  omega  graphs  exhibited  similar 


behavior  to  those  of  0kland  (1970) . 


The  forecasts  of  sea  level  pressure,  temperature,  and  wind  fields 
were  not  significantly  different  in  terms  of  RMS  error.  That  is, 
the  divergent  initialization  did  not  improve  the  forecast  of  these 
variables.  The  SI  forecasts  for  sea  level  pressure  were  best  for  the 
unbalanced  experiments,  especially  for  the  forecast  initialized 
directly  with  observed  data. 

The  precipitation  forecasts  were  compared  in  several  ways.  The 
forecasts  for  Forecasts  3,  A,  and  5  were  scored  with  the  objective 
procedure  presented  in  Appendix  3.  This  procedure  is  a  type  of 
pattern  recognition  evaluation  and  can  determine  if  the  forecast 
precipitation  is  in  the  correct  location  when  compared  with  the 
observations.  When  measured  with  this  technique,  the  forecast  with 
divergent  initial  conditions  did  much  better  in  the  first  three  hours. 
However,  over  the  entire  forecast  period,  all  three  forecasts  scored 
equally  well. 

The  total  precipitation  amounts  over  the  entire  domain  were 
compared  with  the  observed  amounts.  It  was  found  that  the  amounts 
for  the  nondivergent  experiments  gradually  rose  until  about  hour  9 
and  then  gradually  decreased.  Although  a  maximum  occurred  in  the 
divergent  forecast  at  about  hour  9,  the  precipitation  amounts  were 
relatively  uniform  for  the  entire  forecast  period.  The  changes  in 
forecast  amounts  for  the  nondivergent  experiments  were  related  to 
the  adjustment  occurring  in  the  model.  In  the  forecast  with 
divergent  initial  conditions,  the  model  was  much  closer  to  a  balanced 
state  at  the  beginning  of  the  forecast. 


Forecast  5  produced  significantly  more  precipitation  in  the 
12-hour  forecast  period.  Much  of  the  improvement  in  the  forecast 
precipitation  amounts  came  in  the  first  three  hours.  The  nondivergent 
forecasts  started  with  a  precipitation  deficit  and  did  not  make  up 
that  deficit  later  in  the  forecast.  The  best  short-term  precipitation 
forecast  was  obtained  with  the  divergent  initialization. 

The  Initial  forecast  precipitation  rates  were  studied  at  selected 
points  and  for  that  portion  of  the  domain  containing  the  highest 
precipitation  rates.  The  initial  precipitation  rates  were  nearly 
uniform  and  the  area  enclosed  by  specific  contours  was  relatively 
constant.  Although  the  highest  observed  initial  precipitation  rates 
were  underforecast,  both  in  terms  of  maximum  rate  and  the  area  covered 
by  a  given  contour,  the  forecast  model  did  indeed  retain  a  substantial 
portion  of  the  initial  divergence. 

Divergent  initialization  on  the  mesoscale  was  successful  in 
improving  the  precipitation  forecast  for  this  data  set.  It  improved 
the  forecast  in  two  ways:  (1)  the  initial  precipitation  rate  was 
almost  uniform  and  allowed  significantly  higher  precipitation  in  the 
first  three  hours,  and  (2)  the  total  precipitation  amounts  through 
the  entire  forecast  period  were  more  realistic. 


8.0  SUMMARY  AND  CONCLUSIONS 


Numerical  weather  prediction  in  general,  the  initialization  of 
primitive  equation  models,  and  previous  attempts  at  the  initialization 
of  the  divergent  component  of  the  horizontal  wind  were  reviewed  in 
this  thesis.  Most  prior  divergent  initializations  were  on  the 
synoptic  scale  and,  on  that  scale,  had  little  effect  on  the  initial 
precipitation  rate. 

Scale  analyses  were  conducted  to  establish  the  form  of  the 
vertical  velocity  and  divergence  equations  appropriate  on  the  mesoscale. 
The  effect  of  the  assumptions  used  in  neglecting  the  time  dependent 
terms  in  the  vertical  velocity  and  divergence  equations  was  analyzed. 

The  scale  dependence  of  divergent  initialization  on  the  synoptic  and 
mesoscales  was  examined.  For  the  vertical  velocity  equation  on  the 
synoptic  scale,  the  diabatic  term  is  of  the  same  order  as  differential 
vorticity  advection  or  the  Laplacian  of  temperature  advection.  Also, 
the  Coriolis  term  was  important.  On  the  mesoscale,  the  diabatic 
term  dominated  in  precipitation  areas  and  must  be  included.  The 
Coriolis  term  was  not  nearly  as  important  on  the  mesoscale  as  it  was 
on  the  synoptic  scale.  For  the  divergence  equation  on  the  synoptic 
scale,  the  Coriolis  term  was  on  the  same  order  as  the  other  most 
significant  terms.  However,  on  the  mesoscale,  the  effect  of  the 
earth's  rotation  was  less  important. 

The  general  divergent  initialization  scheme  required  the  solution 
of  elliptic  partial  differential  equations  on  a  limited  domain. 
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Since  the  domain  was  not  cyclic  or  periodic,  boundary  conditions  on 
the  dependent  variables  were  required.  After  an  examination  of  what 
properties  the  boundary  conditions  should  possess,  suitable  boundary 
conditions  for  geopotential,  omega,  and  stream  function  were  stated. 

A  new  method  for  the  determination  of  boundary  condition  on  velocity 
potential  on  a  limited  domain,  mathematically  equivalent  to  a  Green's 
function  solution  for  an  infinite  domain,  was  found  such  that  the 
low-wavenumber  boundary  variation  of  velocity  potential  was  accurate. 
This  method's  effectiveness  was  demonstrated  on  analytic  velocity 
potential  fields. 

The  forecast  model  used,  the  PSU  mesoscale  model,  and  its 
nondivergent  initialization  scheme  were  discussed  briefly.  The 
synoptic  case  chosen  for  study  was  also  presented. 

The  finite-difference  form  of  the  omega  equation  was  derived. 
From  experiments  conducted  with  this  equation,  it  was  concluded: 

(1)  The  quasi-geostrophic  omega  equation  with  a  diabatic 
term  overestimated  omega,  especially  in  precipitation 
areas. 

(2)  The  diabatic  term  in  the  mesoscale  omega  equation 
was  the  single  most  important  term. 

(3)  The  largest  single  uncertainty  in  the  determination  of 
omega  and  the  divergent  wind  component  was  the 
representativeness  of  precipitation  observations. 

This  uncertainty  was  transmitted  to  the  divergent  wind 
components  themselves.  The  uncertainty  was  large 


enough  that  there  is  no  advantage  to  calculating 
vertical  velocities  directly  in  sigma  coordinates 
in  order  to  avoid  interpolation  error. 
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The  balance  equation  consistent  with  the  model  was  derived  in 
finite-difference  form  in  sigma  coordinates.  The  procedure  for 
applying  the  balance  equation  was  outlined. 

Five  experiments  were  conducted  which  started  from  the  same 
synoptic  time  but  which  were  initialized  in  a  different  manner. 

The  initializations  were: 

(1)  independent  (unbalanced)  analyses  of  winds  and 
temperatures ; 

(2)  as  in  (1)  except  the  divergence  was  removed  from 
the  winds; 

(3)  nondivergent  winds,  geopotential  calculated  on 
pressure  surfaces  from  the  balance  equation,  and 
hydrostatically  derived  temperatures; 

(4)  nondivergent  winds,  geopotential  calculated  on 
sigma  surfaces  from  the  balance  equation,  and 
hydrostatically  derived  temperatures;  and 

(5)  nondivergent  and  divergent  wind  components, 
geopotential  calculated  on  sigma  surfaces  from  the 
balance  equation,  and  hydrostatically  derived 
temperatures. 
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At  the  completion  of  each  balanced  initialization,  superadiabatic 
lapse  rates  were  present  in  the  lowest  model  layer  in  the  vicinity 
of  the  surface  low  pressure  area.  These  were  removed  by  applying  a 
convective  adjustment  procedure  to  the  balanced  data. 

The  forecasts  were  compared  and  the  following  conclusions  were 
reached : 

(1)  There  was  no  large  difference  in  the  external  gravity- 
wave  noise  characteristics  of  the  forecasts  as  measured 

1^*1  ,3P* 

by  graphs  of  |-r— — ]  and  | — 5-  versus  time.  However,  the 
3t  3t 

forecasts  with  unbalanced  initial  conditions  had  a 
consistently  higher  noise  level  than  the  balanced 
forecasts. 

(2)  There  was  no  significant  difference  in  the  RMS  forecast 
errors  of  sea  level  pressure  or  temperature  and  the  wind 
components  at  the  400-,  500- ,  and  700-mb  levels.  This 
result  was  expected.  The  SI  scores  for  sea  level 
pressure  were  best  for  the  experiments  with  unbalanced 
initial  conditions. 

(3)  For  the  experiments  with  balanced  initial  conditions, 
the  precipitation  forecasts  were  scored  by  an  objective 
procedure.  For  the  first  three  hours  of  the  forecast, 
the  divergent  initialization  experiment  scored  higher 
than  the  other  forecasts.  For  the  entire  12-hour 


period,  all  the  forecasts  scored  equally  well. 


(4)  The  nondivergent  forecasts  experienced  a  precipitation 
deficit  in  the  first  half  of  the  forecast  period 
because  the  model  had  to  develop  a  divergent  component. 
They  produced  nearly  the  same  amount  of  precipitation 

as  the  forecast  with  divergent  initial  conditions  during 
the  last  half  of  the  forecast  period.  Therefore,  the 
divergent  initialization  produced  the  most  accurate 
precipitation  forecast  for  this  case. 

(5)  The  precipitation  amounts  for  each  hour  were  more 
realistic  in  the  forecast  with  divergent  initial 
conditions.  The  nondivergent  forecasts  displayed  an 
oscillation  associated  with  the  mutual  adjustment  of 
the  model's  nondivergent  and  divergent  wind  components. 

(6)  The  model  "remembered"  the  divergent  component  in  the 
experiment  with  divergent  initialization.  The  initial 
precipitation  amounts  were  nearly  uniform.  The  area 
covered  by  the  highest  observed  initial  precipitation 
rates  was  under forecast  and,  therefore,  the  divergent 
initialization  experiment  did  not  predict  as  much 
precipitation  as  observed. 

(7)  Since  the  experiment  with  divergent  initial  conditions 
was  initially  closer  to  a  balanced  state  and  its 
initial  precipitation  rate  was  the  most  realistic,  we 
conclude  that  the  diagnosed  omega  values  were  sufficiently 
accurate  to  be  useful  in  the  divergent  initialization 


performed  here. 


This  thesis  has  the  limitation  that  only  one  data  case  was 
studied.  Even  though  this  divergent  initialization  experiment 
significantly  improved  the  short-term  precipitation  forecast  for  this 
data  set,  it  cannot  be  concluded  that  divergent  initialization  will 
improve  the  precipitation  forecast  in  every  case.  However,  the 
divergent  initialization  procedure  that  was  developed  was  general 
and  should  produce  improved  short-range  precipitation  forecasts  on 
other  cases. 

8.1  Suggestions  for  further  research 

The  balancing  procedure  described  in  Chapter  3  which  is  currently 
used  by  the  PSU  model  should  be  improved.  It  can  be  made  more 
consistent  with  the  model  while  still  balancing  on  pressure  surfaces. 
Presently,  temperatures  are  derived  at  dot  points  and  then  averaged 
to  get  temperatures  at  cross  points  where  temperature  is  defined  in 
the  model.  Temperatures  should  be  derived  directly  at  cross  points 
since  every  unnecessary  averaging  or  interpolation  step  introduces 
error  in  an  initialization. 

Manually  digitized  radar  (MDR)  data  should  be  used  to  obtain 
the  precipitation  rate  required  by  the  divergent  initialization  in  the 
omega  equation  diabatic  term.  These  data  are  described  in  Moore  et  al. 
(1974) .  MDR  data  can  provide  a  precipitation  rate  which  is  more 
spatially  representative  than  raingauge  observations,  since  MDR  data 
represent  entire  grid  squares.  MDR  data  can  also  be  obtained  more 
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readily  than  rainfall  observations.  Unfortunately,  MDR  data  are 
available  only  over  the  eastern  half  of  the  United  States  and  therefore 
are  of  only  limited  usefulness. 

The  divergent  and  nondivergent  initializations  could  be 
improved  by  incorporating  a  three-dimensional  analysis  scheme  before 
the  balancing  procedure.  This  analysis  should  consider  significant 
levels  reported  in  rawinsonde  data.  This  change  should  improve  the 
vertical  consistency  between  the  mass  and  momentum  variables,  eliminate 
most  of  the  superadiabatic  lapse  rates,  and  therefore  enable  the 
model  to  produce  better  forecasts. 

The  precipitation  forecast  produced  by  the  divergent  initializa-  ; 

I 

tion  presented  here  should  be  compared  with  a  dynamic  initialization  < 

on  the  same  data  set.  This  would  determine  which  technique  produces 

the  better  precipitation  forecast.  It  is  possible  that  the  divergent  I 

initialization  would  produce  comparable  forecasts  and  at  a  significantly 

reduced  computer  cost. 

The  divergent  initialization  procedure  should  be  applied  on  the 
synoptic  scale.  Then  the  results  could  be  compared  with  those  of 
other  investigators. 

For  the  forecasting  of  significant  precipitation  events,  better  i 

i 

forecasts  would  result  if  a  complete  radiation  parameterization  were  • 

incorporated  into  the  model.  This  would  be  most  useful  in  convective  j 

situations. 

The  divergent  initialization  scheme  presented  here  should  be 
extended  to  include  the  time-dependent  terms  in  the  omega  and 


divergence  equations.  The  divergent  initialization  scheme  could  then 
be  used  at  grid  increments  as  small  as  10  km.  At  a  10-km  grid 
increment,  an  improved  moisture  initialization  coupled  with  a 
divergent  initialization  would  perhaps  yield  even  better  precipitation 
forecasts. 

The  divergent  initialization  scheme  presented  here  should  be 
tested  on  a  grid  increment  of  about  60  km  for  a  heavy  precipitation 
event.  For  example,  the  scheme  should  be  quite  useful  at  that  grid 
increment  in  the  initialization  of  hurricane  models.  If  precipitation 
rate  data  is  available  (e.g.,  from  satellite  data)  for  a  hurricane 
case,  this  technique  should  provide  a  realistic  hurricane 


initialization. 
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APPENDIX  1 


DERIVATION  OF  THE  FINITE-DIFFERENCE  FORM  OF  THE  OMEGA  EQUATION 
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For  each  term, 
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we  will  now  derive  its  FD  analog, 
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We  will  solve  the  FD  analog  of  (2.23)  for  omega  at  full  (standard- 
pressure)  levels  except  for  the  top  (200  mb)  and  the  bottom  (1000  mb) 


levels.  Therefore,  o  is  interpolated  to  o  using 

Si.j,k'+4s  i.j.k 

o  .  Since  we  do  not  have  temperature  data  at  the  lowest  pressure 
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level,  we  set  o.  ,  .  *  o.  .  Therefore,  T1  becomes,  for  level  k, 
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The  terms  on  the  RHS  of  (2.23)  are  forcing  functions.  Since  we  do  not 
yet  know  omega,  we  do  not  know  or  v^.  Hence,  we  will  first  describe 
the  terms  we  can  compute:  T3,  T5,  T10,  and  Til.  Also,  the  forcing 
functions  will  first  be  computed  at  half-levels. 
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(A1.5) 

+  m  ,  v  (T  ,  .  )  )  ] 

t,i  *1>j>k  i.j.kyyy 


»’ 
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where  Q  is  the  heating  rate  per  unit  mass.  To  determine  the  vertical 

It 

distribution  of  Q,  we  follow  Krishnamurti  (1968a)  and  Anthes  .  The 
specific  humidity,  q,  is  defined  as  the  ratio  of  mass  of  water  vapor 
to  the  total  mass  of  air  volume.  If  q  is  changing  with  time  following 
a  parcel ,  then 


—  «  -  l 

dt  dt 


(A1.6) 


is  the  heating  rate  per  second  per  unit  mass  of  air  (L  -  latent  heat 
of  condensation  580  cal/g  of  water).  We  assume  that,  in  areas  of 
precipitation,  the  air  is  saturated  with  respect  to  water  vapor.  In 
that  case,  expanding  the  material  derivative  on  the  RHS  of  (A1.6) 
yields 


dh 

dt 


3qs 

L<-nr  +  YH*v«s 


3qs 

+  w  sr* 


(A1.7) 


We  further  assume  that,  in  regions  of  precipitation,  the  first  two 
terms  on  the  RHS  of  (A1.7)  are  negligible  compared  to  the  third. 
That  is. 


•ft 

Anthes  (1976)  assumed  a  parabolic  distribution  of  omega  in  the  vertical 
arising  from  the  latent  heating.  The  entire  portion  of  this  description 
dealing  with  the  determination  of  the  parabolic  omega  profile  is  taken 
from  his  unpublished  notes. 
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dh 

dt 


-  Lu 


3p 


(A1.8) 


Integration  of  (A1.8)  over  the  depth  of  the  precipitating  column  per 
unit  horizontal  area  can  be  written 


I  (Pt  dh  .  d2  L  fPt  3qs  . 

II  dFdp-d?-lj  “3Tdp  • 


(A1.9) 


We  can  calculate  4?»  the  heating  rate  through  the  column  per  unit  area, 
at 

2 

from  the  rainfall  rate,  R  (cm/d).  First  convert  R  to  mass/ (cm  s): 


cm  H20  d  1  ■  g  cm  ^  1^0  d  1  x  1  d/8. 64x10^  s 
-  1.1574xl0"5  g  cm-2  s-1  H20 
LR  »  6.71x10  3  cal  cm  2  s  1 


Therefore, 


281  joule  m  2  s  3 


4^  -  LR  -  281  R 
at 


Combining  (A1.9)  and  (A1.10)  yields 


Pt  3qs  ,  g  dQ 

“  w  dP  L  d? 


(A1.10) 


1.134  R  g  s*3  H20  for  [R]  -  cm  d-1  (Al.ll) 


1.134xl0-6  R  cb  s"1 


l 

) 

1 


'  <  **/ 


i  ( 


Pressure 


1 


Wt 


Fig.  35.  Parabolic  convective  omega  profile. 


Now,  we  know  t —  from  the  initial  data  since  we  know  the  temperajture 

dp 

at  the  pressure  surfaces. 

The  third  assumption  in  the  determination  of  Q  is  that  the  omega 
* 

profile  that  results  is  quadratic  in  pressure  (reference  Fig.  35). 


io(p)  »  a  +  bF  +  cp 


The  boundary  conditions  on  the  omega  profile  are 


0)(pg)  -  <ut 


(A1.12) 


(A1.13a) 


ui(p t )  -  0 


(A1.13b) 


, «  V  '■V  ■* , 


We  also  know  that 


w(p  )  “  u> 

'*m  max 


Therefore,  we  can  write 


ai  ■  a  +  bp  +  cp 
t  *s  s 


tii  *  a  +  bp  +  cp 
max  m  m 


0  -  a  +  bpt  +  cpt 


Eliminating  a  between  (A1.14a)  and  (A1.14c)  yields 


2  2 

ut  “  b<ps"pt)  +  c(ps  ~pt  ) 


Eliminating  a  between  (A1.14a)  and  (A1.14b)  yields 


“t  "  “max  “  b(ps-pm>  +  c(pg2-pm2)  ' 


(Al. 13c) 


(Al. 14a) 


(Al.  14b) 


(A1.14c) 


(Al. 15a) 


(Al. 15b) 


Dividing  (Al.l5a)  by  (p  -p.)  and  (A1.15b)  by  (p  -p  )  and  eliminating 

sc  s  m 


b  between  the  resulting  equations  yields 


0)  -U)  0) 

_ t  max  t 

(ps-pm)(pm-pt)  ~  ps~pt 


(Al. 16) 


#jr- 
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We  now  get  b  In  terms  of  c  from  (A1.15b): 


b 


U)  -U) 

t  max 

-  c(p  +p  ) 

p  -p  *3  rm 

*s  m 


(A1.17) 


Given  b  and  c,  we  obtain  a  from  (A1.14c).  Then  (A1.12)  can  be  solved 
for  omega  at  any  pressure  level. 

To  calculate  c,  we  need  to  know  u>  .  We  write  (Al.ll)  in  FD  form 

max 

as 


KMAX 

where  R'  =  1.134x10  b  R  and  KMAX  is  one  less  than  the  number  of  levels. 
Substituting  (A1.12)  into  (A1.18)  yields 

KMAX  KMAX  KMAX  , 

kfx  a6qk4^  +  bpk4^  6V*  +  CP  6qk^  "  R'  ’  (A1‘19> 

Equation  (A1.19)  can  be  rewritten 

af x  +  bf 2  +  cf ^  -  R'  (A1.20) 


where 


KMAX 

fl  “  1  6qk-^  ’  (A1.20a) 

k"l  ™ 


{ 

I 

I 

I 

l 

I 

I 

I 

1 

I 

1 

l 

f 

I 

I 

I 

I 

I 

I 


KMAX 

f2  ’  Pk-Hs  5qk+i5 
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(Al. 20b) 


and 


KMAX 

f3  “  Pk+Js  5qk+Js  *  (Al. 20c) 

Note  that  we  can  calculate  f^,  f ^ >  an<^  ^3*  Solving  (A1.20)  for  c  yields 
R'  -  af  -  bf, 

c  .  - - £  .  (A1.21) 

r3 


Now,  solving  (A1.14c)  for  a  yields 


2 

a  ■  -  bpt  -  cpt  .  (Al. 21a) 

Using  (A1.17)  in  (A1.21a)  for  b  and  inserting  that  result  into  (A1.21) 
leads  to 


w  -u 

c  -f  LR’  +  ([-US 
f3  ps'pm 


-  c(Pg+Pm)]  pt 


cpt  }  fl  *  [ 


U>  -<i> 

t  max 


p  -p  -  c(Ps+pm)]  f2] 


s  m 


(A1.22) 


to  —co 

R'  +  (P,.^  +  *2> 


Let  h  be  the  denominator  of  (A1.22).  Now,  substituting  (A1.16)  into 

Pg-Pm 

(A1.22)  and  solving  for  w  using  the  result  that  -  • 

max  P3-Pt 


•  j 


1 


v 


1 

i 

i 

} 

i 

; 

i 

< 

i 


1/2,  we  get 


*  “t  ‘  a)t(pm'pt)(ptfrf2)  R'  (p3-pm)  (pm'Pt) 

oo  h  h  .  (A1.23) 

max  - - - 

1  "  “IT1  (ptfl+f2) 


To  summarize,  we  can  calculate  from  observations  R',  f^,  f2>  f^, 

pg,  oos>  and  upon  setting  pt>  we  know  p^.  Therefore,  we  can  use  (A1.23) 

to  get  u>  .  We  then  use  (A1.22)  to  get  c,  (A1.17)  to  get  b,  and 
max 

(A1.14a)  to  get  a.  Omega  as  a  function  of  pressure  is  now  given  by 
(A1.12).  Since  we  know  the  omega  profile  due  to  latent  heating,  we 
can  return  to  (Al.ll)  and  calculate 


L  Sp, 


(A1.24) 


where 


ak+4 

Pk+»s 


(A1.25) 


0.622  „  ,,,  AL  r  1  1  ■, 

-  x  0.611  x  exp  L-rrr-  -  - - J 

Pk-Ms  R*  273  Tk+*s 


Finally,  given  Qk+ij»  we  calculate  T10  using 


t(Wxx  +  (Wyy] 

CP  Pk+is 


(A1.26) 


Til:  Til  -  f  3p“ 
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■W* 


I 

I 

I 

I 

1 

I 

I 

I 

! 

i 

I 

f 


First  of  all,  this  term  only  applies  In  the  lowest  level  (K»KMAX  +  1). 
Secondly,  It  only  applies  to  the  lowest  10  cb  of  the  PBL.  Hence  Ap 
in  Til  is  10  cb.  As  in  the  forecast  model, 

Frx  “  8  Ps  CD  VKMAX+1  “KMAX+l  ^A1' 

where  Fr  is  the  x-component  of  the  frictional  force,  p  is  the 
X  s 

-3  -3 

surface  air  density  (here,  10  g  cm  ) ,  =  0.002,  Ug^^^  is  tlie 

1000-mb  u  component,  and  *-s  the  1000-mb  wind  speed,  calculated 

from 


VKMAX+1 


(uKMAX+l  +  VKMAX+1) 


(A1.28) 


Similarly, 


Fr 


y 


o  o  C  V  v 

8  Ms  D  KMAX+1  KMAX+1 


(A1.29) 


Recall  that  we  now  know  the  forcing  functions  T3,  T5,  T10,  and 
Til  at  half-levels.  Therefore,  we  interpolate  them  to  the  full  levels 
from  K  -  2  through  K  =  KMAX.  Let  F  represent  the  total  forcing 

J  >  K 

function  at  each  grid  point.  We  can  now  rewrite  Tl,  T2,  and  F  ,  , 

y  J  >  ^ 

using  (A1.2a)  and  (A1.3),  as 


m,  .[ (a  a),  ,  ,  )  +  (o  uk  .  ,  )  ] 

si,j  k  j  i*j  »k  yy 


(A. 30) 


~  ,  6u> 

+  “i.j  f i, j  [(Vi,j,k)x  "  ^i.j.^y3  6^  5pk  =  Fi,j,k  ' 
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Mathematically,  this  is  a  second-order  elliptic  PDE  and,  given 
boundary  conditions  on  the  top,  bottom,  and  side  faces  of  the  domain, 
(A1.30)  can  be  solved  by  a  standard  relaxation  technique  (Haltiner, 
1971).  The  appropriate  boundary  conditions  that  will  be  used  are 
given  in  Section  2.3.3.  Let  R^  be  the  residual  (the  difference 
between  the  LHS  and  RHS  of  (A1.30))  after  the  ntk  iteration  using 
«“  .  .  .  Then,  we  want  to  find  ,  such  that  ,  *  0.  First, 

leaving  out  map  factors,  we  write 


i.J.k 


(Ax) 


2  [(os“  -“i+l.j.k  +  (os“  ^i-l.j.k 


+  (a  wn) 


s  'i,j+l,k 


+  (a  ti)n) 


a  'i.j-l.k 


“vVj.k1 


+  r“l..i.k-l  “l.J.k 

Pk-l~Pk+l  pk-l"Pk 


(A1.31) 


>k _ i  i  i  » k+l  i  _  p. 

pk~pk+l  i,j,k 


where  B  . 

i.J.k 

gives 


fi,j  Ci,j,k' 


Adjusting  w 


n+1 
i>  j  »k 


such  that  R 


n+1 

i.j.k 


0 


(Ax) 


2  [(os“  )i+l,j,k  +  (CTs“  ^i-l.j.k  +  (asU  ^i, j+l,k 


+  ^s^i.j-l.k  *  4<Vn+1)i,j,k] 


n  n+1  n+1  n 

+  -  LLj.  ~  “i.J.k  _  ~ 

pk-l~pk+l  pk-l"pk  Pk-Pk+1 


(A1.32) 


.  ■  V 


Subtracting  (A1.31)  from  (A1.32)  and  simplifying  yields 


n+1  n 

“i.J.k  ”  “i.J.k 


(A1.33) 


i.j.k 


(— V-  +  ) 

ix)  ^k-1  ^k+1  ^k-1  ^k  ^k+1 


After  applying  (A1.33)  repeatedly,  we  will  have  omega  at  every  grid 
point  in  the  30  x  35  x  7  domain.  With  this  knowledge,  we  can  apply 
the  method  for  the  determination  of  x0  discussed  in  Section  2.3.4  to 

D 

determine  u  and  v  at  the  half  levels.  Interpolation  provides  u  and 
XX  X 

v^  at  the  full  levels.  Now  we  can  return  to  calculate  the  remaining 
forcing  functions  in  (2.23).  They  are  T4,  T6,  T7,  T8,  and  T9. 


T4:  This  term  is  exactly  like  T3  except  in  (A1.4)  replace  u 


u.  ,  v,  ,  and  v, 

^i.j.k+l  ^i, j ,k  ^i, j ,k+l 


i  >  j  >k 

with  u  ,  u  ,  v  , 

xi,j,k  Xi,j,k+1  xi,j,k 


and  v  ,  respectively. 

Xi,j,k+1 


T6:  This  term  is  exactly  like  T5  except  in  (A1.5)  replace  u 


i,  ,  v.  ,  and  v, 

'*'i,j,k+l  *i,j,k  *ij,k+l 


Ki,J.k 


with  u 


Xi> j  >k  Xi,j,k+1  xi,j,k 


and  v 


,  respectively. 


4,j,k+l 


(A1.34) 


D ,  .  ,  ji  is  available  from  the  determination  of  u  and  v  since 


t2  *i.i. 


6“k 

fc+*s  "  Di,j.k+5S  ”  "  6p. 


Therefore,  we  can  write  T7  as 


6? 


f.  .  D 


i,  j  i.I.k+Jj  6p 


1.-1. fr* 


k+h 


T7  is  valid  at  half  levels. 


T8‘  '!?“§ 


T8  «  f  ,  ~~  Z?  -r^~  C .  .  . 
i.j  <5pk  k  6pk  i,j,k 


T8  applies  at  full  levels. 


I9:  f  If  ’“'V  '* 


T9  ■  f .  ,  7^—  [u?  .  .  v  -  w?  u 

i,j  5pk  i,j,kx  5pk  i,j,ky  5pk  *1(jf 


T9  applies  at  full  levels. 


(A1.35) 


(A1.36) 


]  (A1.37) 
k 


To  incorporate  these  terms  into  the  omega  relaxation  scheme,  we 
must  first  interpolate  T4,  T6,  and  T7  to  full  levels.  Note  that  T8 


i 


^  -«i  -T 


and  T9  already  apply  at  full  levels.  Second,  only  one  of  the  above 

terms  contains  «...  itself.  Expanding  (A1.36)  we  get,  for  the  nfct 
1*3  >* 

Iteration, 


rl  /  n 

IT 


+  *  ) 


pk+VPk-%  2  1,j,k+1  1,3  ,k  pk+l“pk 


(A1.38) 


1  .  n  .  n  .  (Ci.1.k  "  5i.i.k-l\ 

2  (wi,3,k +  “i.j.k-^  ’  vCr  ’ 


We  can  now  place  this  term  on  the  LHS  of  (Al.30),  perform  the  appropriate 
modification  to  (A1.31)  and  (A1.32),  and  in  place  of  (A1.33)  we  get 


n+1  n  _ i*1  ,k _ 

Ui,j,k  *  “i.j.k  T1  +  T2  +  T3  +  T4 


(A1.39) 


where 


jLIA 

.  _s  2 


- i-J-lil-  ( - i -  +  - i - ) 

Pk-l~Pk+l  Pk-l"Pk  pk~pk+l 


fi,J(Ci,3,k-l  ~  ‘i.J.k?, 
(pk-l_Pk+l) (pk-l“pk) 


f i , j(Cl.J,k  ~  Ci,J.k+ll 


(pk-l'pk+l) (pk"pk+l) 


•  •* 
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After  applying  (A1.39)  for  a  sufficient  number  of  iterations  over 
the  domain,  we  have  omega  defined  over  the  domain  to  the  accuracy 
permitted  by  (2.23)  and  the  observations. 

Given  omega  over  the  domain,  we  can  compute  the  divergence  and 

velocity  potential,  the  latter  of  which  provides  u  and  v  .  This 

XX 

completes  the  set  of  data  required  for  the  divergent  initialization 
on  sigma  surfaces. 


t 


I 


i 


l 


i 


„  * A  ■*».' 
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APPENDIX  2 

DERIVATION  OF  THE  FINITE-DIFFERENCE  FORM  OF 
THE  BALANCE  EQUATION 

The  term  designators  used  refer  to  the  terms  in  (6.6), 


T1  and  T2: 


i  i 

'T.  .  *  2Ax  (*i+Ss,j-^  +  ^i-Jj.j+%  "  <t’i-+^,i-JS  "  ‘‘‘i-b.j-^ 

**  J 


where  i  and  j  refer  to  dot  points. 


*?■  +  ♦i.j^  ■ 


+  4Ax’  j  (4>i,j-+4s  +  *1-1,  _  *i,j-%  *1-1, J-V 


T1  *  <*x'  }x.  "  [  A4AJA^  (*i+l,j+l  +  ♦i.j+l  '  *l+l,j  '  *i,j> 

1-9  J 


+  ^Ax^  ^i.j+l  +  *1-1, j+1  “  *i,j  “  *i-l.J> 


^  ^  (*i+l,j  +  *i,J  "  *1+1. j-l  '  *l,j~l) 


^Ax^  ^  (*i,j  +  *1-1, j  '  *i, j-1  “  *i-l,J-l)] 


where  i,j  now  refer  to  cross  points. 


Similarly,  for  T2 


-»W' 


J.*/ 


"  ife1  +  ♦i^,J-«S 

J 


where  i,j  refer  to  dot  points. 


T2  *  (<ti  X  )  m—iiA  [  (4  +  4  -A  -*  ) 

'y.  Ax  1  4  Ax  wi+l,j+l  9i+l,j  9i,j+l  9i,J; 

J 


+  —L^hAzh.  (j.  +*  _  4  _  4  ) 

+  4Ax  ^i+l,j  9i+l,J+l  9i,j  9i,J-r 


"  ^Ax^  (*i,j+l  +  *i,j  “  *i-l,j+l  ^i-l.J5 


^Ax1  ^  (*i,j  +  *i,j-l  "  *1-1, J  ~  *i-l,j-l)1 


where  i,j  now  refer  to  cross  points. 

Summing  T1  and  T2  yields  the  LHS  of  (6.6): 


T1  +  T2 


JAl 


2  (Ax) 


2  ^mi+*s, j+*j*i+l, j+l  T  “i-^.j+V'i-l.j+l 


+  m. 


+  ai+*5,j-^i+l,j-l  +  “i-^.j-^i-l.J-l 


+  “1+4, ]-h  +  “i-^.j+^i  +  mi-1s,  j-^^i.  j  ^ 


Note  that  in  FD  notation. 


4  m. 


-xy 


i,J  '  "  “i-Hf.J-H*  +  mi+*a.J-»*  +  “i-Jj.j-Ms  ^  “1-4, i-h 


+  m,  ,  .  iL,  +  m. 


The  RHS  of  (6.6)  (terms  T3  through  T12)  are  the  forcing  functions. 

Once  these  are  evaluated,  they  can  be  used  with  the  expression  for 

the  sum  of  T1  and  T2  to  solve  for  <f> .  , . 

*»  J 


T3:  u  -  ^(^  +  u^  where  i,j  refer  to  dot  points. 

Note  that  p*u  »  up*^.  Let  pi  ^  j  and  ^et 

«"*■>  y 

A.  ,  »  p.  .u.  ..  Then  —  I  *  — —  (A,,,  ,  +  2A  +  A  ,) 

i, j  i. j  i, j  ®  1>:j  j  i+l.J  i.3  1-l.J 


T3a  “  m  'i,j  -  8  [mi^j+Jj  (Ai+l,j+>s  +  2Ai,j+»s  +  Vl.jW 


+  «lf jJs  <Ai+l,J-«s  +  Ai,j-»4  +  i-l.j-V 


(u*  T3a)x  -  ^  ^  <Ai+lf  j+i  +  2Ai,j+l  +  Ai-l,j+l) 


+  miJ  (Ai+l,j  +  ^i.j  +  Vl.j5  (ui,j+l  +  Ui,j> 


(ui,j  +  "i.j-l*  (Ai+l,j  +  2Ai,j  +  Ai-l.j) 


+  <Ai+l,j-l  +  2Ai,j-l  +  ^-l.J-l* 


.  *  <  , 


Let  AI.  , 
i»  J 

*  Ai+l,j+l 

+  2Ai,  j+i  +  Ai- 

A2. 

»  A,  ,  + 

2A  +  A,  ,  , 

i.3 

i+l,j 

i.J  i-l,j 

A3 

*  A . 

+  2A  ,  +  Aj 

i.J 

i+i.j-i 

i,j-l  i 

Simplifying  yields 

4  T3a)x  16Ax  1  m  4J,  ni.  .  “i,J 


3J.  [•  faJti ,  -i.il  ai,  j  a2 , 


i,j+l 


(ui,1  +  “i.i-l* 


i.J-i 


(u*  T3a) 

P±.j 


^i.J1 


44.4  +“44> 


16PUAX  “i.j+l 


AI  +  — 1^+1  ~ — 1»J~1  A2 

Ui  mi>:J  i»  j 


(ui,1  +  Vi-i> 


i.J-l 


— —  T3c 
16Ax  i,j 


*3m] 


where  i,j  still  refers  to  dot  points. 


T3i+*S,j+is  “  "  32^)^  (T3Ci+1-J+1  +  T3CiJ+l  "  T3Ci+l.J  "  T3Ci,J) 


Note  that  l+^.j+ij  is  a  cross  point.  For  example. 


l 


^  ' 


^  V  -H.'  -**'  . 
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T3c 


1+1. J+1  ^i+l.J+2  +  ui+l,j+l)  ^ 


i+1, j+1  P1+i( j+i  mi+2,j+2  i+l,j+l 


+  <Ui+l,j+2  ~  ui+l,j)  A2 


m 


i+1, j+1 


i+1, j+1 


■(ui+l,i-~.U-L+i^.  A3  j 
“i+l,j  1+1'J+l] 


T4:  Note  p  v  *  vp. ^ .  Let  p.  .  *  pT3^  and  B  *  p  v 
*  *  *  i*j 

^  x 

Then  m  li,  j  (Bi,j+l  +  2Bi,j  +  Bi,j-1)  * 

(B 


— y 


T4a  -  - 


I  r_  1 

Q  ^ 


m  8  m 


i+^s.j 


i-^.j+l  +  2Bi+*s,  j  +  Bi-^,j-l) 


— i -  (B<  i,  ix,  +  2®,  i,  ,  +  B,  ,  .  ,) 

“i-^.j  i_Ss,;5+1  i_>S»j  i-^.j-1 


mi^1  ,  +  u,  4> 

y  16Ax 


T4a)„  -  ^  [-l^J - ^ 


mi+1>j  [(Bi+l,j+l  +  2Bi+l,j  +  Bi+l,j-l) 


(B,  +  2B  ,  +  B . 


mljj  i, j+1  i.j  i» j~l 


(ui,j  +  ui-lfj)lmlfJ  (Bi, j+1  +  2Bi,j  +  Bi,j-1) 


(Bi“l> j+1  +  2Bi-l,j  +  Bi-l,j-l)] 


Let  B11J  -  B1+1  j+1  +  2Bi+1>j  + 

B2i,j  *  \j+l  +  2Bi,j  +  Bi,j-1 

B3i,j  ”  Bi-l,j+l  +  2Bi-l,j  +  Bi-l,j-l  ' 
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Simplifying  yields 


,  +U,  i> 


(“m.i  '  Vl,i* 


(7  T4a)  -  B1  + 

1  “1X  *1+1,1  *>3  *1,J 


<U1,1  + 


1-1.1 


83l,l'  • 


T4b  ,  - - X4a) 

M  Pt>J  y 


m  i+J  r(Ultl.J  B1  +  — -+Li  B2 

16Pi,JAX  1  »i+ltj  ^  "ij  ^ 


(ui,1  +  Vl,J> 


i-l.J 


777—  T4c,  . 
16Ax  i,j 


b3i,j! 


Therefore,  T4  at  cross  point  i+^i.j-^S  can  be  written 


T4.,^,  - - (T4c . . ,  ...  +  T4c,  ,  -  T4c . . ,  ,  -  T4Cj  .) 

l'+4S,j+iS  32  (Ax) 2  i+1»j+1  i,J+l  i+1»J 


T5:  T5  is  analogous  to  T3  except  that  in  T3c,  u  is  replaced  with 
v  and,  because  the  outside  differentiation  is  with  respect 
to  y  instead  of  x,  the  final  expression  is 


T5,.,  (T3c.  ..  ...  +  T3c . . .  ,  -  T3c.  ...  -T3c,  ,) 

i-Hi.j-^  32 (Ax) 2  i+1»3+l  i+l.j  i,j+l  i.j 


■  ■ : i v>ks£i> ;*  V ;<*  ‘  *>.V  ‘  i ■  in  1  "T- it" 1  iftf  ^ **  '■ 


T6:  T6  Is  snslogous  Co  T4  except  chat  in  T4c,  u  is  replaced 
with  v  and,  again  because  of  the  outside  differentiation, 
Che  final  expression  is 


Tfe  .  .  _ _ l***! Jjji  (rUc  +  T4c  -  T4e  -  T4c  "> 

^i+is.j-Ms  32(dx)2  i+1*j+1  i+l,J  i,j-l  i.r 


i  6  T*y  — 

T7:  P*u) 


1  6  T*y  _o 
m  6o  0  u 


=5?  5?  (a  "P*  } 


Therefore,  T7 


id 

m  15  (°  U  > 


T*y  i  .  .  .  . 

CTi,j,k+4  “  4  ^i+^.j+^.k+^s  +  0i+ij,j-Ji,k4ii  +  °i-4, j+*i,k+*i  +  ai-*S,  j-H.k-tV 


“i.j.k+^s  2  ^ui,  j,k+l  +  Ui,j  ,k^ 


— xy  , 

•  y  ^  «  •  • 

CTi,j,k+*s  Ui,j,k+Jj  "  8  ^i+is, j+^.k+is  +  0i-Hs, j-^.k+Jj  +  °i-*s, j+^.kH-^ 


+  °Fi-1i,J-!4,k+Ss)(ui,j,lcH  +  “i.j,^ 


V 
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T7a 


I(o4 


+  a. 


i.j,k+*s  8m1(J(ak+1  “  CTk)  i-Ms,j-5s,k+l 


+  a 


+  ai-is,j-4,k+l)(ui>j,k+3/2  +  ui,j,kW 


(ai-^,j4is,k  +  ai+i5,j-vs,k  +  ai-H,5-^,k  +  a±-h,i-h,v? 


^Ui,j,k-Hs  +  “i,k,k-^ 


T7b, 


8(ak+i  “  °k}  t-i,**** 


T7  »  -  _ _  />r7V>  j-  T7Vi 

16Ax(ak+1  -  ofc)  lT,Di+l,J+l,fc*ij  i.j+l.Wfs 


T7bi+l,j,k4As  "  T7bi,j,k4i5) 


T8:  T8  is  analogous  to  T7  except  that  the  u  is  replaced  by  v 

in  T7b.  Also,  the  outside  differentiation  is  with  respect 
y  instead  of  x  and  the  final  expression  is 


T8 


(T7b. 


lrt.J-*  16Ax(ok+1  -  ak)  V1'ui+l,j+l,k4fs  +  T7bi+l,j,k-H5 


T7bi,j+l,k4is  "  T7bi,j,k+is)  * 


T9:  T9a 


mRT*7  p* 


— xy  _t\ 

mp*  (1  + 

°P* 


^iT 


p*  (i  +  pt/°) 


Let  P.  ,  -  T?* 


'Zi  “  4  (Ti+*5,j4ii  +  Ti-^,j“54  +  T1-5s,J-^4  +  Ti-Js,j-Ji) 


pi,j  “  4  (p*i+»s,j+is  +  P*i4Js,j-!s  +  P*i->S,j+is  +  p*i-%,j-%) 


_ y 

p*  “  2Ax  (p*i-^s,j+%  +  p*i-*s,j+is  '  p*u%,j-h  ~  p*i-h,j-h> 


Therefore,  T9a  is 


T9a  ,  mi»j  +  Ti+iS.j-4  *  + 

i,j  2  Ax  (P*^^  +  P*14^,j-Jj  +  p*i-J*,J-tf*  +  P*i-JS,j-l5) 


(1  +  pt/o) 


*  "*  -  T9a,^  , +  T9a  4  ,,  -  T9a_  .  -  T9a^  .) 

i+^j-^  4  (Ax) 2  i+1>3+l  i,j+l  i+l.j  i»  j 


T10:  T10  is  analogous  to  T9  except  that  pA  replaces  p*  y and  the 

y  x 

outside  differentiation  is  with  respect  to  y  instead  of  x. 


P*y  “  2Ax  (p*i-t^s, j+^s  +  P*i4^,j-J5  "  p*i-4,j+ls  '  P*i-4, j-JS) 


t 

i 


T10a  „  ^1*1  +  +  V^.W, 

i,j  2Ax  (P*i444>j4is  +  +  P*i-Si,  j-Hi  +  p*iJStJ-lS) 

(p*i^j.i-^  *  P*l+ij.1-%  ~  ~  p*i-%.  1-%^  i 

(1  +  Pt/o) 


T1°i-^J+i5  *  "  4 (Ax) (T1°ai+]  >1+1  +  T1°ai+l.j  "  T1°ai,j+1  "  T1°ai,j) 


Til:  Tlla 


fp*v  _  iv^7?  m  fv 

— xy  — xy  m 

mp*  mp* 


Tlla  ,  »  — | 
i,j  m 


Therefore, 


(Tllai+i,j+i +  Tllai,j+i  -  Tllai+i,j  -  Tllai,j) 


T12:  T12a 


fp*u  _  £up*  ^  ft, 

— xy  — xy  m 

®P*  mp* 


Therefore, 


2Ax  (T12al+l,j+l  +  T12at+l,j  -  T12ai,j+1  “  T12al,j) 


.  :  \ 
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APPENDIX  3 


i 

OBJECTIVE  PRECIPITATION  SCORING  PROCEDURE 


I 

I 

I 

I 


In  this  appendix,  we  present  the  objective  precipitation  scoring 
procedure  developed  by  Anthes  (unpublished) .  Mesoscale  models  are  now 
producing  forecasts  that  are  finer  in  scale  than  synoptic-model 
precipitation  forecasts.  Hence  a  mesoscale  model  might  forecast  the 
correct  intensity  and  shape  of  a  precipitation  field  but  displace  the 
field  by  some  small  distance.  This  forecast  would  receive  a  poor 
score  with  a  conventional  scoring  method  but  the  forecast  does  contain 
useful  information.  The  skill  score  presented  here  was  designed  to 
detect  forecast  displacement  errors. 

Let  PF  and  PO  be  N  by  M  matrices  of  forecast  and  observed 
precipitation,  respectively.  An  example  matrix  of  the  29  by  34  fields 
used  here  is 

P29,l  P29,34 


P  P 

2,1  *2,2 


P  P 

1,1  1.2 


...  P 


1,34 


(A3.1) 
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I 

I 

I 
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The  skill  of  a  forecast  will  be  computed  for  various  shifts  or 
"lags"  of  the  forecast  and  observed  fields.  A  shift  of  the  forecast 
field  a  distance  kAx  to  the  left  with  respect  to  the  observed  field 
will  be  given  by  a  positive  integer  parameter  k;  similarly,  a  downward 
shift  will  be  given  by  a  positive  1.  When  the  forecast  and  observed 
fields  are  shifted,  a  skill  score  will  be  computed  for  the  overlap 
region  only.  Let  A  be  the  total  number  of  points  in  the  overlap 
region.  Then  A  is  given  by 


A  -  (N  -  8.)(M  -  k) 


(A3. 2) 


Define  the  variance  of  PO  in  the  overlap  region  as 


N-i,  M-k  _ 

Z  Z  (PO  -  por  i>0,  k>0 
i-1  J-1  1,3 


where 


PO 


.  N-i  M  _  , 

r  z  z  (po.  -  po)  n>o,  k<o 

i-l  j-l-k  1,3 

.  N  M-k  _  2 

A  Z  Z  (PO  -  PO)  i<0,  k>0 

i-l-i  j-1  1,3 


(A3. 3) 


-1  N  M  —  2 

‘  Z  Z  (PO.  .  -  PO)  i<0,  k<0 


i-1-1  j-l-k  i,J 


PO  -  A-1  Z  Z  PO 
i  j 


(A3. 4) 
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Similarly,  the  variance  of  forecast  precipitation,  S  ,  is  given  by 

PF 

(A3. 3)  except  that  (P0,  -  P0)^  is  replaced  by  (PF,  -  PF)^.  The 

elements  ^  of  a  matrix  C  are  defined  as 


»->-  M-K  _  _  t>0 

<P0‘.J  -po)<PFi«,i+k-pF)  i>o’ 

H  _  _  ,>0 

1!l  j-U  <P°1,J  -  P°KP!'i«,J+k  -  PF)  k<0 


ASP0SPF  cJt,k 


(A3. 5) 


N  M-k  —  —  *<0 

<poi.j  -M)<ppi«,^-pp)  JS- 

N  M  —  _  *<0 
r  E  (P0  -  P0)(PF  -  PF) 

i-l-A  j-l-k  i,J  i+i,j+k  k<0 


Let  £  and  k  vary  over  some  range.  For  example,  let  i  •  k  ■  2. 
Then  the  matrix  C  is 


-2,-2 


-1,-2 


°-l,l 


(A3. 6) 


The  elements  c^  ^  of  the  matrix  C  are  the  correlation  coefficients 
between  the  observed  and  forecast  precipitation  fields  for  a  given 


displacement  defined  by  i  and  k.  The  maximum  element  of  C  represents 
the  highest  correlation  between  the  observed  and  forecast  patterns. 

The  values  of  1  and  k  represent  the  spatial  offset  necessary  to  achieve 
the  best  score.  Good  forecasts  would  have  a  large  maximum  c.  ,  (1.0  is 

I  K 

perfect)  and,  for  the  maximum  c£  k>  low  values  of  *  and  k  (0  and  0  are 
perfect) . 


■  — 
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